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SECTXCM  I 
INTRODUCTION 


HULL  is  a  system  of  computer  programs  which  solve  the  finite-difference 
analogs  of  the  partial  differential  equations  descriptive  of  a  non-conducting, 
inviscid,  elastic/plastic  continuous  medium  subjected  to  dynamic  loading. 

HULL  is  capable  of  running  cylindrical  or  cartesian  coordinate  calcula¬ 
tions  in  two  spatial  dimensions  or  cartesian  coordinate  calculations  in  three 
dimensions.  HULL  is  a  multimaterial  Eulerian  code  which  employs  a  diffusion 
limiter  to  preseve  material  interfaces.  Up  to  ten  materials  can  be  handled  in 
a  single  calculation.  A  mixture  of  any  of  these  materials  is  permitted  in  any 
one  cell.  Solid  material  properties  are  modeled  by  the  Mie-Gruneisen  equation 
of  state  (Reference  1).  Plasticity  is  modeled  by  the  Von  Mises  yield  criteria 
which  can  be  modified  by  plastic  work  hardening  and  thermal  softening.  Slip 
surfaces  are  automatically  handled  at  material  and  failure  interfaces. 
Failure  can  be  determined  by  simple  stress,  strain,  or  tri-axial  criteria. 
Failure  is  modeled  by  inserting  a  void  or  low  pressure  gas  in  the  failed  zone 
and  allowing  the  tensile  stresses  to  relax.  Recompression  of  the  failed  zone 
is  not  inhibited. 

Large  computational  meshes  can  be  handled  by  HULL  since  only  a  small 
portion  of  the  mesh  need  be  in  central  memory  as  the  calculation  proceeds. 
Time  snap-shots  at  arbitrary  intervals  can  be  retained  on  permanent  disk  or 
tape  files  for  subsequent  analysis  or  problem  restart.  The  HULL  system  has 
several  service  routines  for  keeping  track  of  this  data  storage  and  for  pro¬ 
viding  graphical  presentation  of  most  features  of  interest. 

Provisions  are  available  to  link  HULL  to  the  Lagrangian  TOODY  code 
(Reference  2)  or  EP1C3  code  (Refezence  3)  for  calculations  in  two  and  three 
dimensions,  respectively. 

Development  of  HULL  (Reference  4)  started  in  1971  at  the  Air  Force 
Weapons  Laboratory  (AFWL).  This  initial  version  was  built  to  solve  nuclear 
weapons  effects  problems  with  special  emphasis  on  nuclear  fireball  rise  and 
air  blast  propagation.  In  1975  a  conventional  weapons  effects  version  of  HULL 
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(Reference  5)  wai  developed  at  the  Air  Force  Armament  Laboratory  (AFATL). 
This  later  version  of  the  code  incorporated  new  physics  modules  for  calcula¬ 
ting  hydrodynamic  and  elastic/plastic  forces.  In  addition,  the  code  architec¬ 
ture  was  modified  to  permit  vectorization  on  pipe-line  machines  such  as  the 
GRAY.  HULL  has  been  continuously  developed  by  the  AFATL  and  AFWL  since  1975 
with  primary  management  by  AFATL.  Previous  versions  or  offshoots  of  HULL  are 
not  currently  supported  by  the  AFATL/AFWL  team. 
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SECTION  II 
THE  HULL  SYSTEM 

The  HULL  system  of  codes  is  composed  of: 

(1)  KEEL  -  an  initial  condition  (or  mesh)  generator, 

(2)  HULL  -  the  finite  difference  code, 

(3)  PULL  -  the  graphics  package, 

(4)  BOW  -  the  dump  tape  library, 

(f.  PLANK  -  the  alternate  input  generator  for  SAIL 

(6)  CONTROL  -  an  interactive  program  status 

monitor,  and 

(7)  SAIL  -  the  system  file  manager  and  executive 

processor. 

All  programs  or  codes  in  the  HULL  system  are  maintained  on  a  single  SAIL 
program  library  file.  Compilation  of  any  program  element  of  HULL  entails 
execution  of  SAIL  to  process  the  program  library,  select  the  proper  subrou¬ 
tines  and  options  desired  for  the  particular  application,  and  produce  a  file 
which  will  serve  as  input  to  the  compiler.  This  processing  activity  is  done 
to  keep  storage  requirements  and  central  processor  time  requirements  to  the 
practical  minimum  consistent  with  code  comprehension  and  structure  flexi¬ 
bility.  The  use  of  SAIL  for  -unning  HULL  and  its  associated  service  programs 
will  be  addressed  in  this  document.  Other  documents  which  describe  the  HULL 
code  are  limited  to  specific  features  or  elements  of  HULL  (References  6,  7, 
and  8).  Use  of  SAIL  for  program  alteration  and  file  maintenance  is  covered  by 
Reference  8. 

Figure  1  illustrates  the  interactions  of  the  program  elements  of  the  HULL 
system.  Each  of  the  programs  are  described  more  completely  in  the  respective 
sections  of  this  report. 

The  general  setup  and  running  of  a  HULL  calculation  (assuming  that  an 
executable  SAIL  module  and  the  HULL  system  file  are  available  on  the  computer 
system)  starts  with  the  formulation  of  an  input  deck  for  KEEL.  This  input 
deck  contains  all  the  specifications  for  formulating  the  initial  conditions 
for  the  calculation.  The  KEEL  input  is  first  read  by  PLANK.  PLANK  processes 
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the  deck  to  determine  the  mesh  size  and  options  to  be  used  for  this  run. 
PLANK  then  generates  a  disk  inpat  file  for  SAIL  which  details  the  structure  of 
program  KEEL.  SAIL  uses  this  disk  file  to  process  the  HULL  system  library 
file  and  assemble  the  KEEL  compiler  input.  The  compiler  input  is  compiled  and 
the  resulting  code  is  loaded  and  executed  as  program  KEEL.  KEEL  reads  the 
same  input  deck  used  for  the  generation  of  KF.FI.  gPFI.  then  produces  the  HULL 
problem  initial  conditions  and  puts  the  results  on  logical  unit  4.  Unit  4  is 
either  a  physical  tape  or  a  permanent  file  disk  device  which  will  store  the 
data  for  running  HULL  (or  PULL).  A  second  optional  tape  file  output  is 
produced  if  detailed  time  history  data  at  several  arbitrary  HULL  mesh  points 
are  desired.  This  is  called  the  station  tape,  and  it  resides  on  logical  unit 
9.  The  station  tape  is  required  for  linking  HULL  to  the  Lagrang  ian  TOODY  or 
EPIC3  codes. 

The  output  data  from  KEEL  on  unit  4  is  used  as  input  to  PLANK  to  help 
SAIL  tailor  the  HULL  program  for  a  particular  run.  This  tailoring  process  can 
be  modified  by  the  HULL  input  deck  as  will  be  described  later.  PLANK  and  SAIL 
produce  the  compiler  input.  When  HULL  is  compiled  and  loaded,  it  reads  the 
initial  conditions  from  unit  4  and  starts  the  calculation.  Snap-shots  of  the 
proceeding  calculation  at  pre-determined  times  (from  input)  are  put  on  unit  4 
for  later  analysis  or  restart.  Station  data  is  also  accumulated  and  recorded 
on  unit  9  if  it  is  present.  Unit  4  snap-shots  or  dumps  are  also  made  periodi¬ 
cally  every  60  central  processor  minutes  to  protect  against  unscheduled 
machine  failures.  Thus,  unit  4  and  unit  9  store  a  continuous  record  of  the 
running  calculation  that  is  sufficient  to  restart  the  calculation  at  any  of 
the  available  unit  4  dump  times. 

Unit  4  and  unit  9  are  the  permanent  record  of  the  calculation  for  subse¬ 
quent  analysis  and  editing.  PLANK  uses  these  tape  files  and  a  PULL  input  deck 
to  instruct  SAIL  in  the  generation  of  the  PULL  code.  PULL  produces  an  ex¬ 
haustive  set  of  graphical  outputs  which  can  be  selected  for  problem  analysis. 
The  PULL  code,  after  being  generated  by  SAIL  and  compiled  by  the  operating 
system,  reads  unit  4  or  unit  9  and  produces  graphics  output  for  the  selected 
plotting  devices. 


5 


SECTION  in 

PSOGKAI  SELECTION  AND  CONTROL 


HULL  System  elements  are  selected  and  configured  by  the  SAIL  code  with 
instructions  from  the  input  data  stream  and  PLANE.  The  default  option  direc¬ 
tory  on  the  SAIL  library  file  is  used  to  simplify  this  procedure  by  remem¬ 
bering  the  primary  option  set  values.  The  HULL  system  default  options  are  all 
keyed  from  the  installation  parameter  INST  and  have  the  value  and  installation 
definition  in  the  following  table. 


Value  of 

INST  Installation  Definition 

1  Air  Force  Weapons  Laboratory  (Kirtland  AFB,  NM) 

2  Air  Force  Armament  Laboratory  (Eglin  AFB,  FL) 

3  McDonnel  Douglas  Astronautics  (Huntington  Bch,  CA) 

4  Army  Ballistic  Missile  Defense  Agency  (Huntsville,  AL) 

5  Science  Applications  Incorporated  (Huntsville,  AL) 

6  Royal  Armament  Research  and  Development  Establishment 
(Fort  Halstead,  Kent,  UK) 


7  Army  Ballistic  Research  Laboratory  (Aberdeen,  MD) 

8  Vought  (Dallas,  T3) 

9  OTI  (Shalimar,  FL) 


Once  the  installation  parameter  is  set,  the  options  which  choose  the  compnter 
system  specific  parameters  are  selected  automatically. 

After  the  correct  options  have  been  assembled  to  tailor  the  HULL  programs 
for  the  installation  machine  and  operating  system,  we  can  further  control  SAIL 
processing  to  produce  the  specific  form  of  KEEL,  HULL,  or  PULL  needed  to  run  a 
given  problem.  The  dump  tape  (unit  4)  produced  by  KEEL  is  the  major  control 
for  setting  up  programs  HULL  and  PULL.  Most  of  the  input  to  KEEL  sots  these 
control  parameters  through  a  unit  4  record  called  the  ZBLOCK.  The  names  of 
the  ZBLOCK  parameters,  the  meaning  of  their  various  values,  end  the  defaults 
are  given  below. 


ZB LOCK  Name 

Values 

Definition 

(Alternate  name) 

(Default) 

PROBLEM  (PROB) 

Any  number  of  the  form 

Arbitrary  problem  identifier 

mi.YYYY 

AREF 

True.,  (False.) 

Y(j=l)  boundary  reflective  (3D) 

ATMOS 

Integer  1 

tropical  atmosphere 

2 

temperate  atmosphere 

3 

arctic  atmosphere 

4 

exponential  atmosphere 

(5) 

constant  pressure  and  density 

BREF 

.TRUE.,  (.FALSE.) 

bottom  boundary  reflective 

Y(j=l),  2D; 

Z(k=l) ,  3D 

BURN 

Integer  (0) 

no  burn  routine 

1 

burn  front  advanced  by  fluid 
state 

2 

programmed  burn 

CODE 

Integer  (1) 

normal  hydrodynamic  -  elastic/ 
plastic  formulation 

2 

include  interactive  dust 

COLD 

(.TRUE.),  .FALSE. 

radiation  cooling  if  false 

CYCLE 

Real 

number  of  time  steps  from 
initial  conditions 

DIMEN 

(2),  3 

number  of  spatial  dimensions 

DT 

Real 

time  step  (seconds) 

ELC 

Real 

energy  last  check,  current 
total  mesh  energy 
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EOS 


1 


Doan-Nickel  air  equation  of  state 


2 

(6) 

10 

ETH  Real 

EXPAND  Real  number  <1(.10) 

FAIL  Integer  (0) 

1 

FLDXER  Integer  0 

1 

(3) 

FREF  .TRUE.,  (.FALSE.) 

GEOM  1 

(2) 

HOB  Real  (0) 

IMAX  Real  (100) 

IQ  Real  (IMAX-1) 

ISLAND  ~  (0) 

1 

2 


Constant  gamma  equation  of  state 

Multi-material  Mie— Gruneisen 
equation  of  state 

tabulated  air  equation  of  state 

theoretical  total  me~H  energy 

mesh  expansion  fraction  for 
rezone 

no  failure  model 

simple  stress  criteria  (SIRAIN=0) 
simple  stress  and  strain  criteria 
(STRAIN  =  1) 
no  fluxing 

single  material  fluxing 

volume  fluxing  with  species  mass, 
volume,  and  energy  retained 

Y ( j=JMAX )  boundary  reflective 
(3D) 

cartesian  coordinates 

cylindrical  coordinates  (2D  only) 

weapon  burst  height  in  cm 

number  of  zones  in  x  coordinate 
direction 

number  of  active  x  coordinate 
zones 

no  effect 

internal  mesh  reflective 
boundary 

reflective  boundary  with  rigid 
body  translation  (not  implemented 
in  3D) 
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JMAX 

Real  (200) 

number  of  zones  in  y  coordinate 
direction 

JQ 

Real  (J MAX-1) 

number  of  active  y  coordinate 
zones 

KMAX 

Real  (1) 

number  of  zones  in  z  coordinate 
direction 

KQ 

Real  (KMAX-1) 

number  of  active  z  coordinate 
zones  (3D) 

LREF 

(.TRDE.),. FALSE. 

X ( i=IMAX )  boundary  reflective 

MAGFLD 

(0) 

no  magnetic  field 

1 

dipole  field  and  MHD  formulation 

METHOD 

0 

no  calculation  in  Lagrange  step; 
used  for  testing  fluxer 

(2) 

(second  order  in  time) 

MLC 

Real 

mass  last  check,  current 
total  mass  in  the  mesh 

NH 

Integer 

number  of  variables  per  zone 

(5  +  3*  NM  +  NHIST; 

2D) 

calculated  by  KEEL 

(6  +  3*  NM  +  NHIST; 

3D) 

NHIC 

Integer  (set  by  PLANK) 

number  of  words  in  core  to 
retain  mesh  variables 

NHIST 

Integer  (set  by  PLANK) 

number  of  fluxed  material  de¬ 
pendent  histories  per  cell 
(stress  components,  plastic 
»crk,  etc.) 

NM 

Integer  (1) 

number  of  materials 

NVARST 

NOP 

Integer  (IS,  2D; 

22.  3D) 

Integer  (0) 

number  of  variables  per  station 

number  of  trace  particles  and 
stations 

NPLPB 

Integer  (calculated 
by  PLANK) 

number  of  planes  per  record 
block  on  TAPE4  (3D) 
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t  . 


r 

NPP 

Integer  (2,  2D; 

3,  3-D) 

i' 

number  of  parameters  per  particle 
(usually  coordinates) 

NSOVPB 

Integer  (calculated 
by  PLANE) 

number  of  rows  per  record  block 
on  unit  4 

NSTN 

Integer  (0) 

number  of  stations 

PTSTOP 

Real  (1.E20) 

problem  stop  time 

RAD 

Integer  (0) 

no  radiation  routines 

1 

equilibrium  radiation  diffusion 

: 

2 

non-equilibrium  radiation 
diffusion 

RADLOSS 

Real 

total  energy  radiated 

l  REZONE 

(0) 

no  rezone 

1 

1 

Horizontal  and  vertical  -  shock 
triggered 

1 

2 

particle  triggered  -  horizontal 
and  vertical 

3 

horizontal  shock  follower 

4 

vertical  shock  follower 

5 

vertical  particle  follower 

6 

same  as  S  with  shocks  diffused 

7 

continuous  (only  rezone  for  3-D) 

8 

vertical  material  2  follower 

RREF 

.TRUE.,  (.FALSE.) 

X(i=IMAX)  boundary  reflective 

STRAIN 

Integer  (0) 

no  strain  components  saved 

? 

1 

strains  saved  (only  if  STRESS  >0) 

I  STRESS 

0 

no  stress  (elastic/plastic 
formulation  not  included) 

j 

(1) 

stress  calculation  included 

1  SOME 

Real 

total  energy  lost  by  cooling 

l 
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TE8AD 

Real 

total  energy  loat  by  radiation 

TLC 

Real 

time  of  laat  check  on  energy  and 
mass 

TREF 

.TRUE. ,  (.FALSE.) 

Y ( j  = J MAX ) ,  2D;  Z(k-KMAl) ,  3D 
boundary  reflective 

TTIME 

Real 

total  problem  central  processor 
time  (seconds) 

ITS  TOP 

Real  (100.. 

total  central  processor  time 
stop  (hours) 

UREZ 

Real  (10.) 

z  coordinate  velocity  rezone 
trigger  (cm/sec) 

vise 

Integer  (0) 

no  explicit  artificial  viscosity 

1 

explicit  artificial  viscosity 
included 

VREZ 

Real  (10.) 

y  coordinate  velocity  rezone 
trigger  (cm/sec) 

VOIDS 

Integer  (0) 

no  voids  allowed 

1 

implicit  voids  pexmitted 

2 

explicit  voids  whenever  material 
VOID  is  found 

WORK 

Integer  (0) 

no  plastic  work  hardening 

1 

plastic  work  hardening  included 

WREZ 

Real  (10.) 

z  coordinate  velocity  rezone 
trigger  (3D)  (cm/ sec) 

XOB 

Real  (0.) 

x  coordinate  of  burst  (3D) 
(centimeters) 

YOB 

Real  (0.) 

y  coordinate  of  burst  (3D) 
(centimeters) 

YGND 

Real  (-1.E20) 

location  of  ground  for  rezone 
down  control;  BREF  is  set  to 
zero  when  mesh  bottom  reaches 
this  value 

11 


The  ZB LOCI  variables  are  identified  on  TAPE4  by  their  names  in  alphabetic 
form  along  with  their  corresponding  values.  KEEL  will  default  ZBLOCK  varia¬ 
bles  to  the  values  indicated  above.  If  a  different  value  is  desired,  i*  "an 
be  input  through  the  KEEL  or  HULL  input  data  stream. 
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SECTION  IY 
PROGKAM  KEEL 


Program  KEEL  defines  the  initial  ZBLOCK,  the  mesh  coordinates,  and  the 
material  property  specifications  (P,  n,  v,  w,  I,  M,  V)  on  a  cell  by  cell  basis 
for  the  entire  mesh.  Inpnt  to  KEEL  is  free  format  with  variable  names  and 
values  delimited  by  blank,  comma,  and  equal  (=).  The  KEEL  input  must  start 
with  the  word  KEEL.  This  allows  program  KEEL  to  search  for  and  find  its  data 
independent  of  its  location  in  the  input  file.  The  KEEL  input  is  logically 
divided  into  three  parts:  initial  ZBLOCK  definition,  mesh  generation,  and 
material  generation.  The  logic  flow  of  KEEL  is  indicated  by  Figure  2. 

A.  ZBLOCK  DEFINITION 

The  ZBLOCK  is  literally  defaulted  for  every  variable  name.  The  only 
variables  that  must  be  set  are  to  define  the  number  of  materials  and  their 
names.  It  is  usual  to  associate  a  problem  number  to  identify  the  type  of 
problem  being  run,  a  user  numeric  identifier,  and  to  use  the  numerals  to  the 
right  of  the  decimal  point  to  indicate  the  sequence  or  series  of  the  calcula¬ 
tion  when  running  parameter  studies.  If  the  BOW  program  library  is  being  used 
to  keep  track  of  problems  and  their  associated  tapes,  it  is  imperative  that 
all  problem  numbers  be  unique.  HULL  problems  can  be  assigned  problem  numbers 
from  00000.0001  to  99999.9999. 

After  the  problem  number  is  defined,  it  is  necessary  to  determine  the 
number  of  zones  to  be  used  in  each  coordinate  direction.  The  HULL  code  is  not 
sensitive  to  moderate  variations  in  cell  size  from  zone-to-zone  in  any  given 
coordinate  direction  (10  percent).  However,  it  is  not  tolerant  of  large 
aspect  ratios  (greater  than  2.5-3).  The  number  of  zones  should  be  chosen  to 
provide  enough  resolution  to  properly  define  material  interfaces.  Choosing 
the  best  number  of  zones  in  each  coordinate  is  always  a  compromise  most 
heavily  weighted  by  economic  considerations  since  the  computer  time  required 
to  run  a  calculation  is  inversely  proportional  to  the  size  of  the  smallest 
zone  and  directly  proportional  to  the  total  number  of  zones  in  the  calcula¬ 
tion.  Therefore,  decreasing  the  zone  size  by  a  factor  of  2  in  a  3D  calcula- 

l 
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START 


Figure  2.  Logic  Flow  in  Program  KEEL 


tion  would  result  in  a  factor  of  16  increase  in  the  run  time.  The  number  of 
zones  in  the  x,  y,  and  z  coordinate  directions  is  set  by  IMAX,  JMAX,  and  KJtAX 
respectively.  Two-dimensional  calculations  will  not  use  KMAX  (the  z  coordi¬ 
nate)  . 

The  next  order  of  business  is  to  select  the  number  of  materials  and  their 
names.  Table  B-l  (Appendix  B)  lists  the  materials  currently  recognized  by  the 
HULL  system.  Any  material  listed  in  Table  B-l  may  be  input  to  KEEL,  but 
present  code  structure  limits  their  total  number  in  a  given  calculation  to  ten 
(i.e.,  NM^IO).  Each  material  defined  must  be  assigned  a  number  from  1  to  NM. 
The  order  is  arbitrary. 

If  deviatoric  stresses  are  not  to  be  included,  set  STRESS  =  0.  Otherwise, 
the  calculation  will  assume  all  materials  have  strength.  If  it  is  desired 
also  to  retain  the  strain  components  for  each  cell,  set  STRAIN  -  1.  This  is 
necessary  if  STRAIN  or  P/T  criteria  are  to  be  used  for  failure.  Similarly,  if 
material  failure  is  to  be  explicitly  modelled,  it  is  necessary  to  set  FAIL  = 
1.  This  has  effect  only  for  STRESS  =  1.  If  FAIL  =  1,  STRESS  =  1.  and  STRAIN 
=  0;  failure  will  be  determined  by  a  default  stress  criteria;  maximum  princi¬ 
pal  >  2.5  yield  stress.  If  STRAIN  =  1,  both  the  stress  criteria  and  a 
principal  strain  >  .7  will  initiate  failure.  These  are  default  values  and 
should  be  used  with  caution  because  of  our  current  poor  understanding  of 
failure  initiation. 

. 

Rezone  choice  may  be  made  by  indicating  the  number  from  0  to  7  for  the 
particular  rezone  selected.  REZONE  =  0  implies  no  rezone.  This  csn  be  al¬ 
tered  later  as  HULL  is  rur  ing.  Only  REZONE  =  7  is  supported  in  three  dimen¬ 
sions.  This  rezone  (in  2D  or  3D)  can,  in  general,  cause  a  continuous  transla¬ 
tion  of  the  mesh  coordinates  in  any  direction.  It  is  currently  implemented  to 
translate  the  entire  mesh  at  a  continuous  velocity  in  the  vertical  direction 
(y  direction  in  2D;  z  direction  in  3D).  The  translation  velocity  is  set  by 
the  zone  which  contains  station  or  particle  number  one.  Use  of  the  continuous 
rezone  may  reduce  the  numerical  method's  implicit  viscosity  to  an  unacceptable 
level.  In  any  case  it  is  best  to  run  most  dynamic  solid  calculations  with 
explicit  artificial  viscosity  by  setting  VISC  =  1. 
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If  it  is  desired  to  run  with  Lagrangian  tracer  particles,  the  variable  NOP 
should  be  set  to  a  value  greater  than  the  number  of  particles  to  be  generated. 
The  particle  routine  will  count  and  set  NOP  to  the  actual  number  inserted.  If 
stations  are  to  be  generated,  NSTN  should  be  set  to  a  number  greater  than  the 
number  to  be  inserted.  Since  all  stations  are  also  particles,  NOP  must  be 
greater  than  the  number  of  particles  plus  the  number  of  stations.  If  a 
physical  tape  is  to  be  used  to  store  station  data  (and  neither  BOW  nor  the 
permanent  files  are  being  used),  then  it  is  necessary  to  set  STAPE  =  VSN, 
where  VSN  is  the  volume  serial  number  of  the  tape.  This  will  only  function  on 
CDC  machinery  and  results  in  the  automatic  issue  of  the  tape  request  by  KEEL, 
Otherwise,  you  must  specify  the  station  tape  unit  9  by  means  of  control  cards. 

The  stability  factor  STABF  (fraction  of  the  Courant  condition  time  step) 
is  currently  defaulted  to  ,7j.  The  default  cycle  one  time  step  is  1  x  10~^ 
seconds.  If  the  problem  type  being  run  has  had  difficulties  starting  during 
previous  HULL  runs,  DT  should  be  reduced. 

For  most  circumstances  an  initial  time  step  DT  of  10-10  seconds  should 
suffice.  If  STABF  has  to  be  set  to  much  less  than  0.5,  something  is  seriously 
wrong;  check  your  initial  conditions. 

With  the  current  version  of  HULL  (Version  117),  ISLAND  =  2  and  REZONE  =  1, 
2,  3,  4,  5,  6  have  not  been  implemented  in  the  three-dimensional  HULL  code. 

If  the  problem  being  generated  has  a  mesh  region  that  is  inactive  (no 
velocity,  high  pressures,  or  explosive  burn  region),  some  initial  problem 
computer  time  can  be  saved  by  setting  IQ,  JQ,  and  KQ  (3D).  This  will  cause 
the  mesh  in  that  part  of  the  coordinate  beyond  X(IQ),  Y(JQ),  or  Z(KQ)  to  be 
left  out  of  the  calculation  until  the  activity  is  such  that  the  velocity  in 
any  of  the  active  region  limit  exceeds  UREZ,  VREZ,  or  WREZ  for  the  X,  Y,  and  Z 
coordinate  directions,  respectively.  This  feature  can  not  be  employed  in 
conjunction  with  the  continuous  rezone  (if  REZONE  =  7).  As  the  activity 
reaches  the  limits  of  the  active  region,  the  corresponding  activity  index  IQ, 
JQ,  or  KQ  will  be  incremented  one  zone  per  cycle  until  it  reaches  the  mesh 
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limit  IMAX-1,  JMAX-1,  or  KMAX-1.  The  defaults  for  IQ.  JQ.  and  IQ  are  IMAX-1, 
J MAX-1,  and  K MAX-1. 

The  last  input  in  this  section  of  KEEL  is  the  header  title.  This  is  80 
columns  of  whatever  information  the  user  wishes  to  use  to  identify  the  calcu¬ 
lation.  This  title  will  be  printed  at  the  beginning  of  the  output  for  each 
HULL  run  and  will  appear  in  the  titles  of  all  plots  produced  by  PULL.  The  80 
columns  of  information  must  be  preceded  by  a  card  image  with  the  word  HEADER. 

The  end  of  ZBLOCK  data  definition  is  signalled  by  the  word  MESH. 

An  example  of  this  part  of  the  KEEL  input  for  a  three-material  problem  to 
calculate  the  stress  induced  in  an  RHA  plate  impacted  by  copper  and  ignoring 
failure  could  be  set  up  by: 

KEEL  PROB=1001 .001 

NM=3  AIR=1  KHA=2  CU=3  DIMEN=2 

STRESS=1  VISC=1  DT=1.E10 

REZ0NE=0  IMAX=80  JMAX=120 

HEADER 

TEST  CALCULATION  TO  ILLUSTRATE  KEEL  INPUT 
MESH 

If  it  is  desired  to  insert  stations  to  record  stress  time  histories  at  10 
points,  then 

NOP=ll ,  NSTN=11 

could  be  placed  between  or  on  any  card  preceding  HEADER.  The  almost  explicit 
(with  aid  of  the  ZBLOCK)  definition  makes  this  portion  of  KF.F.I.  input  relative¬ 
ly  easy  to  use. 

B.  MESH  GENERATION 

The  next  logical  element  of  the  KEEL  input  is  definition  of  the  mesh 
coordinates.  This  part  of  the  KEEL  input  begins  with  the  keyword  MESH. 
Several  options  are  possible.  The  simplest  is  a  linear  mesh  in  all  coordi- 
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nate  directions.  This  is  set  up  by  establishing  the  mesh  boundaries  in  each 
coordinate  direction.  The  variable  names  are: 


XO  coordinate  of  X(I=0) 

YO  coordinate  of  Y(J=0) 

ZO  coordinate  of  Z(K=0)  (3D  only) 
XMAX  coordinate  of  X(I=IMAX) 

YMAX  coordinate  of  Y(J=JMAX) 

ZMAX  coordinate  of  Z(K=KMAX)  (3D  only) 


The  mesh  coordinates  are  established  for  the  boundary  of  each  zone  having 
indices  (1,  J,  K)  at  the  position  midway  between  its  center  and  zones  (I+l,  J, 
K),  (I,  J+l,  K),  and  (I,  J,  K+l).  Thus,  the  definition  of  XO,  YO,  and  ZO  is 
to  locate  left,  back,  and  bottom  of  the  first  cell.  Each  of  these  variables 
must  be  paired  with  a  numeric  value.  There  are  restrictions.  In  two  dimen¬ 
sions  o.  v  X0=0  is  allowed.  In  all  cases  XO,  YO,  and  ZO  are  defaulted  to 
zero.  The  zones  will  be  of  constant  size  for  this  mesh  specification  with 
zom  size: 

DX=(XMAX-XO)/IMAX 
DY=( YMAX-YO) /JMAX 
DZ= ( ZMAX-ZO ) / ZMAX  (3D  only) 

Thus,  for  the  values  of  IMAX=80  and  JMAX=120  in  the  preceeding  sections,  the 
input 


MESH  XMAX =40  Y0=30  YMAX=90 


would  produce  a  mesh  with  DX=DY=0.5  centimeter. 


There  are  usually  economic  constraints  to  running  calculations  with  zones 
small  enough  for  sufficient  resolution  of  all  effects  everywhere  in  the  mesh. 
The  best  compromise  frequently  is  to  concentrate  a  fine  constant  grid  in  the 
region  of  most  interest  and  increase  the  zones  in  all  directions  away  from 
that  region.  This  is  especially  important  for  reasonable  treatment  of  solids 
since  perfectly  transmissive  mesh  boundaries  for  the  HULL  code  have  not  yet 
been  formulated.  Thus,  if  the  calculational  mesh  is  too  small,  disturbances 
will  partially  reflect  from  its  boundary  and  propagate  to  the  interior  of  the 
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mesh.  These  reflected  disturbances  may  disrupt  the  validity  of  the  calcula¬ 
tion  in  the  region  of  interest.  This  mesh  generation  option  is  obtained  with 
the  keywords: 

MESH 

CONSTANT  SUBGRID  NX=NN  X0=X  XMAX=XX 
NY=MM  Y0=Y  YKAX=YY 
NZ=LL  ZO=Z  ZMAX=ZZ  (3D) 

NX,  NY,  and  NZ  are  the  number  of  zones  in  the  constant  subgrid  region  with 
boundaries  XO,  XMAX;  YO,  YMAX;  and  ZO,  ZMAX  respectively.  The  zone  dimensions 
within  this  region  are  defined  by 

DX=(XMAX-XO) /NX 
DY=(YMAX-YO)/NY 
DZ= ( ZMAX-ZO ) / NZ  (3D). 

NX,  NY,  and  NZ  are  defaulted  to  IMAX,  JMAX,  and  KMAX,  respectively.  The  re¬ 
maining  IMAX-NX,  .TMAX-NY,  and  KMAX-NZ  zones  outside  the  constant  region  will 
be  increased  in  size  by  a  constant  factor  as  the  zone  index  increases  beyond 
NX,  NY,  or  NZ.  The  factors  are  RXPOS,  RXNEG;  RYPOS,  RYNEG;  and  RZPOS,  RZNEG 
for  zones  with  indices  greater  or  less  than  the  constant  subgrid  boundary  in 
the  X,  Y,  and  Z  coordinate  directions,  respectively.  The  default  value  is  1.1 
for  all  these  factors,  i.e.,  a  10-percent  growth  rate  will  be  applied  for 
adjoining  zones  as  they  get  progressively  further  from  the  constant  subgrid 
region.  This  is  a  good  value  to  use  in  most  situations.  If  it  is  desired  to 
have  the  exterior  limits  of  the  mesh  stop  near  some  limit,  then  the  variables: 

X0LIM=X  XMLIM=XX 
YOLIM=Y  YHLIM=YY 
ZOLIM=Z  ZMLIM=ZZ 

can  be  specified  along  with  a  large  value  of  the  corresponding  multiplying 
factors.  The  result  in  the  X  coordinate  direction  would  be  to  place  X(0) 
XOLIM  and  X(1MAX)  <.  XMLIM.  The  other  coordinates  would  be  treated  in  a 
similar  manner  for  their  corresponding  variables. 
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If  the  above  lacks  flexibility  for  a  particular  application,  the  user  can 
specify  the  size  of  each  zone  with  the  variable  name/value  pairs 

MESH  NX=il  DX=xl  NX  =  i2  DX  =  x2 

NY=jl  DY=yl  NY=j2  DY  =  y2 

NZ=kl  DZ=zl  NZ=k2  DZ  =  z2 

where  il,  i2  are  the  number  of  zones  with  zone  size  xl  for  the  first  il 
zones,  and  zone  size  x2  for  the  next  i2  zones  until  the  sequence  sum  is  equal 
to  IMAX.  XO  defines  the  origin.  The  mesh  coordinates  Y  and  Z  will  be  handled 
in  an  analogous  way.  The  individual  value  of  each  occurrence  of  NX,  NY,  and 
NZ  can  be  from  one  to  IMAX,  JMAX,  and  KMAX,  respectively.  Thus,  for  a  brute 
force  mesh  setup  one  could  specify 


MESH  Y0=0 

X0=0 

NX=1 

DX=.l 

NX  =2 

DX=  .2 

NX =3  DX=.3 

NY=2 

DY=.l 

NY=4 

DY=.5 

NY=2  DY=1 

would  produce  mesh  coordinates  (assuming  IMAX=6  and  JMAX=8)  of: 

Xj=0,  .1,  .3,  .5,  .8,  1.1,  1.4;  i=integers  [0,  6] 

Yj^O,  .1,  .2,  .7,  1.2,  1.7,  2.2,  3.2,  4.2;  j=integers  [0,  8] . 

The  search  for  additional  mesh  data  will  terminate  with  the  keyword 
GENERATE. 

C.  MATERIAL  GENERATION 

The  KEEL  input  data  starting  with  the  word  GENERATE  is  used  for  inserting 
materials,  particles,  and  stations  in  their  desired  mesh  coordinates.  Al¬ 
though  the  nuclear  options  for  isothermal  spheres  and  other  constructs  still 
remain  in  the  code,  they  have  not  been  used  in  several  years  in  this  version 
of  the  code  and  probably  do  not  execute  properly.  Their  use  is  covered  in 
previous  documents  (Reference  4)  and  will  not  be  repeated  here.  Their  func¬ 
tions  can  be  replicated  by  the  syntax  which  follows. 
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Material  can  be  inserted  into  the  mesh  by  KEEL  from  two  basic  sources — 
from  a  previous  HULL  or  EPIC3  calculation,  or  from  the  !1EE1 .  input  data  stream. 
To  input  data  from  a  previous  calculation  the  keyword  FIREIN  must  be  used. 
This  is  immediately  followed  by  the  word  HULL  or  EPIC3  depending  on  the  data 
source.  The  next  variable  name/value  pair  to  occur  is 

TAPE=VSN  or  PFN=permanent  file  name 

Impermanent  file  identifier 

where  VSN  is  the  volume  serial  number  of  the  tape  which  contains  the  HULL  data 
and  PFN  and  ID  are  the  permanent  file  name  and  file  identifier  (ID),  respec¬ 
tively,  of  the  EPIC3  data.  HULL  data  is  supplied  by  any  previously  run  HULL 
calculation  in  a  standard  TAPE4  format.  EPIC3  data  is  provided  by  a  special 
dump  routine  which  is  described  in  Reference  30.  The  KP.RT.  input  at  this  point 
would  appear  as 

GENERATE 

FIREIN  HULL  TAPE=VSN 


FIREIN  EPIC3  PFN=pfn  ID=id. 


The  construct  for  defining  the  material  and  location  destination  in  the 


mesh  is: 


PACKAGE  material  geom  DELETE  geom  1  geom  2 


I* 

f] 

I; 


where  material  is  the  HULL  name  of  one  of  the  materials  defined  in  the  first 
part  of  the  KEEL  input  data  stream.  This  material  must  also  have  been  present 
in  the  calculation  which  produced  the  HULL  or  EPIC3  data  tape  being  input  or 
this  package  will  have  no  effect.  Geom  is  the  geometric  figure  and  its 
specifications  which  define  the  region  in  the  mesh  that  is  to  be  filled  from 
the  input  data  tape.  The  geometric  figures  which  follow  (geom  1,  geom  2  to  a 
total  of  9)  are  deleted  from  the  first  geometric  configuration.  The  word 
DELETE  is  not  necessary,  since  the  appearance  of  a  second  and  subsequent 
geometry  specification  within  a  PACKAGE  implies  deletion  of  that  region  encom¬ 
passed  by  its  boundaries. 
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Each  of  the  geometric  figures  is  accompanied  by  a  list  of  parameters  which 
define  their  descriptions  and  position  in  the  mesh.  The  following  constructs 
are  for  describing  two-dimensional  material  generation. 

RECTANGLE  Xl=xl  X2=x2  Yl=yl  Y2=y2 

A  rectangle  with  vertex  coordinates 
(xl,  yl),  (xl,  y2) ,  (x2,  yl) ,  (x2,  y2) . 

TRIANGLE  Xl=xl  X2=x2  X3=x3 

Yl=yl  Y2=y2  Y3=y3 

A  triangle  with  vertex  coordinates 
(xl.yl),  (x2,y2),  (x3,y3). 

CIRCLE  XC=xl  YC=yl  RADIDS=r 

A  circle  of  radius  r  with  center  at  (xl,yl). 

ELLIPSE  A=yl  B=b  C=xl  D=d 

An  ellipse  with  center  (xl,yl)  and  semiaxis  length  b 
in  the  y  coordinate  and  d  in  the  x  coordinate 

PARABOLA  A=a  B=b  C=c 

All  (xj.yj);  y j  <  a+b  (x^-c)*. 

HYPERBOLA  A=a  B=b  Oc  D=d 

All  (xj.y j) ;  yA  2  a+b  l+(xi-c)*/d* 

GNRLFIT  A=a  B=b  C=c  D=d  E=e  F=f 

All  (ij.yj);  yj  2  axj*  +  bxj4  +  cxj*  +  dxA*  +  exj  +  f. 

CORVE  TABLE  NPT  =  n  (xl.yl)  (x2,y2) - (xn.yn) 

For  n  <  100  points  and  xl  <  x2  <  x3  <  -  <xn 

^m+l_ym 

All  (xi,yi);  xffl  2  >  *m+i  and  yA  2  - ) 

Tn+l  *m 

The  three-dimensional  geometry  definitions  are  described  in  a  similar  way. 

Box  Xl=xl  X2=x2  Yl=yl  Y2=y2  Zl=xl  Z2=z2 

The  inside  of  the  prism  with  corners  (xl,yl,zl), 

(xl,yl,z2),  (xl,y2,zl),  (xl,y2,z2),  (x2,yl,zl), 

(x2,yl,z2),  (x2,y2,zl),  (x2,y2,z2). 


WEDGE  Xl=xl  Yl=yl  X2=x2  Y2=y2  Y3=y3  X3=x3  Zl=zl  Z2=z2 

The  inside  of  the  triangnlar  prism  with  base  in  the  XY 
plane  and  corners  (xl.yl),  (x2,y2),  (x3,y3)  in  the  planes 
Z=zl  and  Z=z2. 

PYRAMID  Xl=xl  Yl=yl  Zl=zl  X2=x2  Y2=y2  Z2=z2 

X3=x3  Y3=y3  Z3=z3  X4=x4  Y4=y4  Z4=z4 

The  inside  of  the  tetrahedron  defined  by  the  triples 
(xj.yj.zj);  i=l,2,3,4. 

CYLINDER  XC=xl  YC=vl  ZB=zl  ZT=z2  RADIUS=r 

The  interior  of  the  right  circular  cylinder  of 
radius  r  with  axis  parallel  to  the  Z  axis  and  centered 
in  the  XY  plane  at  (xl.yl).  The  cylinder  ends  are  at 
Z=zl  and  Z=z2 . 

ELLIPCYL  A=a  B=b  C=c  D=d  ZB=zl  ZT=z2 

The  interior  of  the  elliptical  cylinder  described  by 
the  parameters  a,  b,  c,  d  as  in  the  two-dimensional 
ellipse  with  ends  at  Z=zl,  Z=z2. 

PARCYL  A=a  B=b  C=c  ZB=zl  ZT=z2 

The  interior  of  the  parabolic  cylinder  analogous  to 
the  ellipse  above. 

HYPERCYL  A=a  B=b  C=c  D=d  ZB=zl  ZT=z2 

The  interior  of  the  hyperbolic  cylinder  analogous  to 
the  ellipse  above. 

SPHERE  XC=xl  YOyl  ZC=zl  RADIDS=r 

The  interior  of  the  sphere  of  radius  r  at  (xl,yl,zl). 

RECTAROT  Xl=xl  X2=x2  Zl=zl  Z2=z2 

The  region  swept  by  the  rectangle  specified  rotated 
about  the  z  axis. 

TRIROT  Xl=xl  Zl=zl  X2=x2  Z2=z2  X3=x3  Z3=z3 

The  region  swept  by  the  specified  triangle  rotated 
about  the  Z  axis. 

CIRCLEROT  XC=xl  ZC=zl  RADIDS=r 

The  region  swept  by  the  specified  circle  rotated  about 
the  Z  axis. 
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ELLJPSEROT  AZ=a  BZ=b  CZ=c  DZ=d 


The  region  swept  by  the  specified  ellipse  rotated 
abont  the  Z  axis. 


PARABLROT  AZ=a  BZ=b  CZ=c 


The  region  swept  by  the  specified  parabola  rotated 
about  the  Z  axis. 


GNRLFIT  A=a  B=b  Oc  D=d  E=e  F=f 


The  region  swept  by  the  polynomial  analogous  to  the 
two-dimens ‘ onal  general  fit  rotated  about  the  Z  axis. 


CURVE  TABLE  NPT=N  xl,  zl.  x2 .  z2 - xn.  zn 


For  n£100,  the  region  swept  by  the  analog  of  the 
two-dimensional  curve  table  rotated  about  the  Z  axis. 


All  of  these  geometric  figure  specification  parameter  names  have  synonyms. 
The  following  names  in  each  set  are  equivalent: 


SET  1 

XO,  Xl, 

XL,  XLEFT,  XC,  XCENT,  XCN1R 

SET  2 

YO,  Yl. 

YL.  YB.  YBOT,  2B0TT0M,  YC, 

SET  3 

ZO,  Zl, 

ZL.  ZB,  ZBOT,  ZBOTTOM,  ZC,  : 

SET  4 

X2,  XR, 

XRIGHT 

SET  5 

Y2 ,  YR, 

YF,  YT.  YTOP 

SET  6 

22,  ZR, 

ZT.  ZTOP. 

An  additional  degree  of  flexibility  is  provided  by  allowing  each  of  the 
above  figures  to  be  displaced  or  rotated  about  each  of  the  coordinate  axis. 
Ihe  displacement  distance  is  specified  by 


XCOxl  YCOyl  ZCC=zl  (3D) 


where  the  occurence  of  one,  two,  or  all  of  the  above  would  displace  the 
geometric  frame  of  reference  by  xl  in  the  X  coordinate  direction,  yl  in  the  Y 
coordinate  direction,  or  zl  in  the  Z  coordinate  direction. 


Potation  of  the  geometric  figure  can  be  accomplished  b / 


AHGLA=a  DANA=da  NDA-~1 

ANGLB=b  DANB=db  NDB=m  (3-D  only) 

ANGLO c  DANOdc  NDOn  (3-d  only) 


Rotation  will  be  accomplished  with  respect  to  the  geometric  frame  of  reference 
defined  by  XCC,  YCC,  and  ZCC.  A  positive  rotation  is  counterclockwise  when 
viewed  along  the  axis  of  rotation.  All  of  these,  any  pair,  or  any  singly 
defined  angle  could  be  defined  for  rotating  about  the  Z,  Y,  or  X  axis  by 
angles  a,  b,  and  c  degrees,  respectively.  The  DANA=da  NDA=1  portion,  when 
present,  indicates  a  replication  of  the  Figure  1  times  at  increments  of  da 
degrees.  All  of  these  variables,  for  specification  of  displacement  or  rota¬ 
tion,  default  to  zero.  An  illustration  using  a  HULL  input  tape  is  shown 
below. 


GENERATE  FIREIN  HULL  TAPE=AA1 
PACKAGE  FE  RECTANGLE  XL=0  XR=10 

YB=0  YT=1 

YCC=-.5  XCC=~5  ANGLA=45 

would  result  in  producing  a  mask  (for  a  2-D  mesh)  in  a  slab,  of  thickness  1 
cm,  centered  on  the  origin,  at  an  angle  of  45  degrees  to  the  X  coordinate 
axis.  Data  would  be  obtained  from  the  HULL  problem  on  TAPE  AA1  and  any  iron 
(FE)  properties  available,  as  seen  through  the  mask,  would  be  inserted  in  the 
new  mesh.  Subsequent  PACKAGE  input  geometries  would  obtain  data  from  the 
same  tape  until  terminated  by  ENDFXRE. 

The  most  common  usage  of  the  package  contruct  is: 

PACKAGE  material  [U=u,  V=v,  W=w,  I=i,  RHO=p] 
geometry  DELETE  geom  1  geom  2  -  geom  9. 

The  definition  of  the  material  and  geometry  parameters  is  as  before.  The 

parameters  enclosed  by  [ - ]  and  the  word  DELETE  are  optional.  Any  of  them 

may  appear.  They  are  intended  to  specify: 

u=X  component  of  material  velocity  (cm/sec) 
v=Y  component  of  material  velocity  (cm/sec) 
w=Z  component  of  material  velocity  (cm/ sec) 
i=internal  specific  energy  (ergs/gm) 
p=material  density  (gm/cm*). 


25 


The  velocity  components  default  to  zero.  The  energy  and  density  default 
to  ambient  conditions  for  standard  temperature  and  pressure,  as  defined  by  the 
equation  of  state  for  that  material. 

Figure  3  is  a  two-dimensiofaal  example  with  a  plot  illustrating  the 
effects  of  the  rotation. 

A  final  use  of  the  geometry  specifications  is  the  generation  of  Lagrangian 
trace  particles.  This  construct  is: 

PARTICLES  Geometry  DELETE  g eon  1,  geom  2, - geom  9. 

The  description  geometry  and  geom  1  -  geom  9  are  treated  as  before,  but  are 
used  to  insert  one  Lagrangian  particle  in  each  cell  whose  center  lies  within 
the  region  defined.  The  word  DELETE  is  optional.  The  particles  are  counted 
as  they  are  inserted  to  keep  the  ZBLOCK  variable  NOP  properly  set. 

The  final  KEEL  input  set  is  station  generation.  This  is  of  the  form 

STATION  XS=xl  x2  x3 - xl 

YS=yl  y2  y3 - ym 

ZS=zl  z2  z3 - zn 

any  number  of  such  sets  can  appear.  They  have  the  effect  of  producing  data 
collection  points  at  each  location  specified  by  xi,  yj,  zk  in  the  list  above. 
This  is  an  implied  loop  that  is  executed  as 

LOOP  z=zl ,  z2 ,  z3  -  zn 

LOOP  y=yl ,  y2,  y3  -  ym 

LOOP  x=xl,  x2,  x3  -  xl 

ENDLOOP  DEFINE  station  with  coordinates  (x,y,z). 

The  maximum  of  1,  m,  n  must  be  less  than  NSTN  as  input  to  KEEL.  The  use  of 
XS,  YS,  ZS  (3-D)  indicates  stations  that  are  to  remain  fixed  with  respect  to 
the  Eulerian  reference.  Lagrangian  stations,  stations  that  follow  the  flow, 
can  be  inserted  by  replacing  the  S  in  XS,  YS,  and  ZS  by  an  L.  A  station  can 
have  Lagrangian  or  Eulerian  behavior  in  each  of  the  coordinate  directions, 
i.e.,  XL,  YS,  ZL  is  permissible. 
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GENERATE 

PACKAGE  FE  CIRCLE  XC=24  YC=1.2  RADIUS=1.2 

RECTANGLE  YB=1.2 

PACKAGE  FE  RECTANGLE  XL=22.8  XR=25.2 

YB=1 . 2  YT=40 

PACKAGE  RHA  V=1 . 5E5 

RECTANGLE  YB=-8.4  YT=0 
XCC=24  YCC=-2.6 
ANGLA=60 


Figure  3.  Plot  Illustration  Rotation 
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After  all  particles,  stations,  and  packages  have  been  processed,  KEEL 
writes  ont  the  mesh  data  on  unit  4  and  the  station  data  on  nnit  9.  The 
formats  are  shown  in  the  accompanying  Figures  4  and  5.  KEEL  also  performs 
checks  for  empty  zones  and  indicates  their  location  and  number.  A  column  list 
and  material  and  particle/station  print  plot  is  also  produced  to  allow  the 
user  to  check  the  results  of  generation. 
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UNIT  4  FORMAT 


IJfV1  iWf}  *?*'«•»  ».V  w»w 


■A. 


(DUMP  TAPE) 


DUMP 


T  555 . ,  PR  OB,  CYCLE,  T  (Header  Synchronization) 

ZBLK(100,2)  (ZBLOCK  Description) 

X(IMAX) ,  YO,  Y(JMAX),  ZO,  Z(KMAX)  (Coordinate  Description) 

H  (NVARPB) -|  (Cell  by  Cell  State  Variables) 

• 

.  NBLKS 

J 

H(NOP*NPP)  IF  NOP>0  and  NOP*NPP  <  NVARPB  (Particle  and  Station 

Coordinates) 

lH(INT(NVARPB/NPP)*NPP)  NOP*NPP  2  NVARPB 


.  Additional  dumps 

• 

666.,  666.,  666.,  666.  (End  of  the  World  Trailer  ID) 

EOF 

EOF 


Figure  4.  HULL  Mesh  Dump  Tape  Format 
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UNIT  9  FORMAT 
(STATION  TAPE) 


f  STAT,  PEOB 

Header  ZBLK  (100,2) 

Informa¬ 
tion  H(N0P*NPP)  NOP*NPP<1020 
otherwise 

LH(NPP*INT(1020/NPP))  for  number  of  records  needed 
fTime  (last  IS  bits  holds  number  of  words  at  this  time) 


HYDRO 

STRESS(2-D) 

STRESS  (3-D) 

PtSTA-last  15 

Density[STA-last  15  bits] 

Density [STA-last 

bits] 

X [Material  Code] 

15  bits] 

Density [Material 

X [Material  Code] 

Code] 

D 

D 

U 

V 

V 

V 

w 

Repeated 

W 

t 

6 

for  next 

active 

if 

V 

station 

till  1024 

* 

word 

SRR-P 

SXX-P 

buffer  is 

SRZ 

SYY-P 

filled 

SZZ-P 

SZZ-P 

SHOOP-P 

SXY 

ERRL 

SXZ 

EZZ 

SYZ 

EH OOP 

EXX 

Y 

EYY 

I(NVARST>15) 

EZZ 

EXY 

EXZ 

EYZ 

Y 

l 

Z 

MATERIAL  CODE  -  bit  number  indicates  material  present  i.e.,  last  15  bits  for 
indicating  materials  bit  on  indicates  presence  of  that  ma¬ 
terial  =  material  3  and  5  is  00000  00000  10100. 

Figure  5.  HULL  Station  Tape  Format 
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SECTION  V 
PROGRAM  HULL 


HULL  solves  the  finite  difference  equations  on  a  mesh  defined  by  KEEL  or 
a  previous  HULL  run.  The  mandatory  input  is  the  initial  conditions  contained 
on  unit  4.  If  station  data  is  being  retained  (NSTN>0),  then  a  unit  9  data 
will  also  be  input.  All  other  input  data  for  running  HULL  will  be  obtained 
from  the  input  data  stream.  The  HULL  input  must  start  with  the  keyword  HULL. 
All  data  on  the  HULL  input  record  is  in  free  format,  delimited  by  blank, 
comma,  and  equal  (=).  The  HULL  input  record  is  logically  divided  into  three 
parts:  restart  information,  ZBLOCK  and  program  control,  and  solid  material 

specifications.  Figures  6  ind  7  are  schematic  representations  of  the  data 
flow  in  program  HULL. 

A.  RESTART  DATA 

The  first  portion  of  the  input  record  immediately  following  the  keyword 
HULL  is  used  to  assure  the  correct  unit  4  input  is  being  read.  The  parameter 
name/value  pair  PROB  =  XXXX.YYYY  is  compared  with  the  header  record  on  unit  4. 
If  they  are  not  the  same  to  within  lO-^,  HULL  terminates  abnormally.  If  this 
test  is  passed,  HULL  checks  to  see  if  the  problem  time  and  cycle  number  on 
unit  4  is  less  than  the  value  desired  for  restart.  If  they  are  not,  the  unit 
4  data  are  read  until  the  condition  is  met  or  the  end  of  file,  as  designated 
by  the  HULL  trailer,  is  reached.  If  end  of  file  is  reached,  HULL  backspaces 
unit  4  to  the  beginning  of  that  dump  or  time  snap-shot  for  restart.  The  time 
and  cycle  for  restart  are  defaulted  to  l.xlOi0.  If  restart  at  a  particular 
time  or  cycle  is  desired,  the  parameter  name/value  pairs:' 

T=t  CYCLE=c 

can  be  inserted  in  the  input  stream  immediately  after  the  PROB=XXXX.YYYY  pair. 
Only  T  or  CYCLE  need  be  input.  The  value  t  instructs  HULL  to  restart  on  the 
first  problem  dump  whose  problem  time  is  greater  than  or  equal  to  t.  Cycle 
starts  are  determined  in  the  same  way  except  the  comparison  is  done  between 
the  value  c  and  '  *LE  as  it  appears  on  unit  4.  To  be  used  for  restart,  all  of 
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ROW/PLANE 

N 

ROW/PLANE 

N-l 

ROW/PLANE 

N-2 

ROW/PLANE 

N-3 

FIRST  WATCH 
SECOND  WATCH 
THIRD  WATCH 
FOURTH  WATCH 


PRINCIPAL  SUBROUTINE  LINKAGE 


FIRST  WATCH 


SECOND  WATCH 


THIRD  WATCH  FOURTH  WATCH 


*— EOSSET 


LstATE 
-AIR 
—  EARTH 
—FIRE 
-WATER 


— BURNHE 
— PTRAN 


L STATION 


‘-HYDRO 


Lstate 

-AIR 

-EARTH 

-  FIRE 

-  WATER 


L-FLUX 


L-FLUX  2 


Figure  7.  Memory  Resident  Logical  Data  Mapping  and  Principal 
Subroutine  Linkage  During  In-Cycle  Processing 
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the  names  PROS,  T,  and  CYCLE  must  appear  before  the  end  of  this  part  of  the 
HULL  input  record.  This  part  is  terminated  by  the  keyword  INPUT. 

A  typical  restart  segment  of  a  HULL  input  deck  is: 

HULL  PR0B=1 .001  CYCLE=105 
INPUT 

This  would  cause  HULL  to  read  the  unit  4  header,  assuring  that  the  tape  con¬ 
tained  information  for  problem  number  1.001  and  advance  unit  4  to  the  first 
dump  whose  cycle  is  greater  than  or  equal  to  105.  If  cycle  100  is  the  last 
cycle  on  unit  4,  it  would  be  positioned  to  restart  at.  cycle  100. 

B.  ZBLOCX  AND  PS 06 KAN  CONTROL 

The  part  of  the  HULL  input  data  following  INPUT  and  preceding  the  end  of 
file  or  the  keyword  MATPRQP  is  used  to  redefine  some  ZBLOCK  parameters  and  to 
specify  program  edit  control  and  tape  dumps.  Any  parameter  value  in  the 
ZBLOCK  can  be  modified  by  this  portion  of  the  input  stream  by  using  the 
appropriate  7BLOCK  name/value  pair  in  free  format  as: 

NAME=value 

where  NAME  is  the  ZBLOCK  name  and  value  is  an  integer,  floating  point,  or 
logical  input  definition,  as  defined  in  program  KEEI-.  Care  should  be  exer¬ 
cised  in  this  usage  since  it  is  possible  to  redefine  the  problem  number, 
current  time,  or  cycle  by  inadvertently  placing  the  keyword  PROB,  T,  or  CYCLE 
after  the  keyword  input  rather  than  before.  Thus: 

HULL  PR0B=1 .001 
CYCLB=24 
INPUT 
T=1 . E-6 

would  result  in  starting  the  calculation  3t  the  first  cycle  2  24  and  then 
redefining  the  problem  time  for  that  cycle  as  l.E-6  seconds.  The  following 
ZBLOCK  parameters  should  not  be  modified  since  they  will  change  the  process 
flow  of  the  code  or  the  meaning  of  an  edited  output  variable. 

34 


Dangerous  HULL  INPUT  ZBLCCK  Redefinitions 


SET  BY  KEEL  TO  SET  BY  HULL  TO 

DEFINE  THE  RUN  PROVIDE  INFORMATION 


CODE 

DIMEN 

I  MAX 

JMAX 

KMAX 

NH 

NHIST 

NM 

NOP 

NROWPB 

NPLPB 

NPP 

NSTN 

NVARST 

PROB 


CYCLE 

ELC 

ETH 

MLC 

MTH 

SUME 

T 

TLC 

TTIME 

TTIME6 

TTIME7 


Resetting  the  values  of  IQ,  JQ,  KQ,  BURN,  FLUXER,  FAIL,  EOS,  or  STABF  will 
produce  results  which  may  be  harmful  to  the  health  of  the  calculation.  They 
should  also  be  used  with  care. 


There  are  several  variables  which  HULL  uses  for  control  during  a  calcula¬ 
tion  that  do  not  appear  in  the  ZBLOCK.  These  are  listed  along  with  their 
function  and  default  values  in  the  tabulation  that  follows. 


Non- ZBLOCK  HULL  Inpmt  Control  Parameters 


NAME  FUNCTION  DEFAULT  VALUE 

DCYCST=n  stop  calculation  in  n  cycles  l.xlO10 

RTSTOP=t  stop  calculation  after  this  run  l.xlO* 

accumulates  t  CPU  hours 


DEBUG=L  print  and  dump  every  cycle  if  L  is  .TRUE.  .FALSE. 

TlMES=n  n=l  produces  36  dumps  per  decade  in  time  1 

n=2  produces  dumps  at  10ms,  30ms; 
every  .Is  from  ,1s  to  1.2s;  every 
0.5s  from  1.5s  to  6.0s 

n=3  produces  dumps  at  interval  DMPINT 
(commonly  used  for  most  conventional 
weapons  problems) 
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DMPINl^At 


e  interval  At  for  producing  dumps 


l.xlO*0 


VMIN=v  minimum  velocity  component  permitted;  1  cm/sec 

if  the  absolute  value  of  any  calculated 
velocity  component  is  less  than  v  it  is 
zeroed 

ITRACE=i  cell  (i,j,k)  will  be  traced  by  printing  l.xlO4 

JTRACE=j  all  intermediate  values  when  SAIL  option 

KTRACE=k  DE3UG=2  is  set 

HEADER  80  columns  from  the  following  card  will 

be  used  to  replace  the  ZBLOCR  header  title 

MRELER  maximum  relative  error  allowed  per  cycle  l.xlO-** 

in  mass  or  energy 

This  section  of  HULL  input  terminates  with  the  keyword  MATPROP  or  an  end  of 
record. 

C.  MATERIAL  SPECIFICATIONS 

All  material  properties  for  HULL  are  defaulted  to  the  best  current  guess. 
If  a  calculation  involves  a  particular  material  specimen  for  which  experimen¬ 
tal  results  are  available,  it  is,  in  general,  best  to  replace  the  default 
values  of  the  yield  strength  parameters  by  their  experimentally  determined 
values.  The  parameters  that  can  be  changed  through  the  material  properties 
specification  input  that  follow  MATPROP  are  in  a  material  dependent  format 
that  can  contain  the  following  optional  name/value  pairs. 

MATPROP 

MAT=n,  YLDST=yl ,  YLDMAX=y2 ,  EPLASI^e,  IFRACl=ifl,  YFRACl=yfl , 
IFRAC2=if2,  YFRAC2=yf2 

ENDPROP 

If  the  keyword  MATPROP  appears,  then  ENDPROP  must  appear.  The  specification 
of  each  set  of  specifications  for  a  given  material  must  be  preceded  by  MAT=n 
where  n  is  the  material  number  to  be  changed.  The  other  variable  name/value 
pairs  are  optional,  but  names  mus£  be  spelled  correctly  or  the  routine  will 
search  forever  trying  to  recognize  them.  In  order  for  EPLAST=e  to  have  an 
effect,  the  ZBLOCK  must  contain  WORK  >1;  otherwise,  the  yield  surface  will  be 
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defined  by  YLDMAX.  When  WORK  is  on  (2D*  plastic  work  (Cp)  is  accumulated  for 
each  zone  which  is  occupied  by  a  solid.  The  yield  surface  is  then  defined  by 
Y  as  a  function  of  the  accumulated  plastic  strain,  as  indicated  in  the  equa¬ 
tion  of  state  description.  Yield  is  alv.ys  a  function  of  internal  specific 
energy.  Thermal  softening  is  accomplished  by  a  tri linear  fit  as  previously 
discussed  in  the  equation  of  state  section.  The  thermal  factor  is  multiplied 
times  the  current  yield  surface  definition  (including  plastic  hardening)  to 
obtain  a  final  value  of  Y  (flow  stress)  against  which  the  yield  test  is  made. 
If  you  do  not  input  these  values,  the  default  values  from  Table  B-2,  Appendix 
B,  will  be  used. 
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SECTION  VI 
PROG IAN  POLL 


Program  POLL  accepts  the  calculational  mesh  data  from  unit  4  or  the  sta¬ 
tion  data  from  unit  9  and  produces  plotted  output  on  a  variety  of  graphics 
display  devices.  The  specific  plots  desired  can  be  specified  by  the  user 
through  the  POLL  data  input  stream.  INPUT  to  PULL  is  free  format  with  input 
parameter  naae/value  pairs  and  keywords  delimited  by  blank,  comma,  or  equal 
(=).  The  PULL  input  must  start  with  the  word  PULL.  Actually  two  different 
program  sets  can  be  produced  by  SAIL  with  the  keyword  PULL.  The  first  is  the 
general  mesh  data  dump  plotter  which  produces  plots  from  the  HULL  unit  4.  The 
second  program  set  is  for  sorting  and  plotting  the  station  tape  data  written 
by  HULL  on  unit  9.  Program  PULL  data  flow  is  depicted  by  Figure  8. 

A.  DUMP  TAPE  PLOTS 

PULL  is  capable  of  making  three  basic  types  of  plots.  A  field  variable 
(density,  pressure,  energy,  etc.)  versus  one  of  the  spatial  coordinates  (x,y, 
or  z)  is  called  a  histogram  and  is  abbreviated  for  use  in  the  keyword  struc¬ 
ture  as  HST.  Lines  of  constant  value  of  a  field  variable  in  a  two-dimensional 
representation  i’e  plotted  for  up  to  20  arbitrary  constants  or  contours. 
These  are  called  contour  plots  and  are  abbreviated  in  the  keyword  calls  by 
CONT.  Plots  of  vector  quantities  with  an  arrow  indicating  vector  direction 
and  the  arrow-length  indicating  vector  magnitude  can  be  plotted  for  all  vector 
quantities  such  as  velocity  and  gradients  of  scalars.  These  plots  are  avail¬ 
able  only  for  two-dimensional  planes.  Vector  plots  are  abbreviated  in  key¬ 
words  as  VECT.  When  gradients  of  scalars  are  selected,  both  a  vector  repre¬ 
sentation  of  the  gradient  and  contour  plots  of  the  gradient  magnitude  are 
produced.  Gradient  plots  are  abbreviated  in  keywords  as  GRAD. 

Histograms  can  be  obtained  in  any  coordinate  direction.  Those  in  the  x 
coordinate  direction  are  called  horizontal  histograms  and  are  abbreviated  by  H 
which  precedes  the  designation  HST.  Histograms  in  the  y  or  z  coordinate 
direction  are  called  vertical  histograms.  They  are  abbreviated  by  V  for 
construction  of  keywords.  Thus,  to  obtain  histogram  plots  of  any  field  vari- 
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Figure  8.  Program  PULL  Data  Flow 


39 


IPgWWPWTOPBWWg  PWPgpB?!  W 


T4-^l,Hll|My- 


able  XXX  the  keyword  would  be  XXXVHST  or  XXXBBST.  Contour  plots  are  called  by 
XXXCONT  and  gradient  plots  are  selected  by  XXXGRAD.  The  allowable  designa¬ 
tors,  XXX  and  their  meanings  are  tabulated  in  Table  1. 


TABLE  1.  BLOT  DESIGNATORS 


POSSIBLE  PLOTS 

XXX _ DEFINITION _ BBST  VBST  CONT  GRAD 


P 

D 

V 

W(3D) 

E 

D 

T 

M 

K 

RP 

RE 

RD 

SXX (3D) 
SYY(3D) 
SRR(2D) 

szz 

SBB(2D) 

SRZ(2D) 

SXY(3D) 

SXZ(3D) 

SYZ(3D) 

S2I 

SP1 

SP2 

EXX (3D) 
EYY(3D) 
ERR (2D) 
EZZ 

EBB(2D) 

ERZ(2D) 

EXY(3D) 

EXZ(3D) 

SYZ(3D) 

E2I 

EP1 

EP2 


HYDROSTATIC  PRESSURE  X 
X  COMPONENT  OF  VELOCITY  X 
Y  COMPONENT  OF  VELOCITY  X 
Z  COMPONENT  OF  VELOCITY  X 
SPECIFIC  INTERNAL  ENERGY  X 
DENSITY  X 
TEMPERATURE  X 
INDIVIDUAL  MATERIAL  DENSITY  X 
SPECIFIC  KINETIC  ENERGY  X 
PRESSURE  IN  ATMOSPHERES  X 


ENERGY  IN  UNITS  CF  AMBIENT  X 
DENSITY  RELATIVE  TO  AMBIENT  X 
X  COORDINATE  NORMAL  STRESS  X 
DEVIATOR 

Y  COORDINATE  NORMAL  STRESS 


DEVIATOR 

RADIAL  STRESS  DEVIATOR  X 
Z  OR  AXIAL  STRESS  DEVIATOR  X 
HOOP  DEVIATORIC  STRESS  X 
SHEAR  STRESS  X 
SHEAR  STRESS  X 
SHEAR  STRESS  X 
SHEAR  STRESS  X 
SECOND  STRESS  INVARIANT  X 
MAXIMUM  PRINCIPAL  STRESS  X 
MINIMUM  PRINCIPAL  STRESS  X 
STRAIN  FOR  SXX  X 
STRAIN  FOR  SYY  X 
STRAIN  FOR  SRR  X 
STRAIN  FOR  SZZ  X 
STRAIN  FOR  SHB  X 
STRAIN  FOR  SRZ  X 
STRAIN  FOR  SXY  X 
STRAIN  FOR  SXZ  X 
STRAIN  FOR  SYZ  X 
SECOND  STRAIN  INVARIANT  X 
MAXIMUM  PRINCIPAL  STRAIN  X 
MINIMUM  PRINCIPAL  STRAIN 


XXX 

XXX 

XXX 

XXX 

XXX 

XXX 

XXX 

XXX 

XXX 

XXX 

XXX 

XXX 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


X 

X 

X 


X 

X 


40 


Thus,  PCONT  would  be  the  allowable  keyword  for  selecting  pressure  contours. 

As  many  of  the  above  as  have  meaning  (i.e.,  SRR  is  not  available  for 
STRESS=0  run)  can  be  entered  for  a  single  FULL  run. 

Contour  plots  can  be  modified  by  a  list  of  variables  following  the  keyword 
to  specify  contour  values  to  be  plotted.  For  example: 

DCONT  =11. 5  2345 

instructs  PULL  to  draw  contour  lines  for  each  of  the  density  values  1,  1.5,  2, 
3,4,  and  5  gm/cm1.  Contours  can  also  be  modified  by  ICTRS=n  and  CONYV=xl , 
x2,  ...,  xn  being  selected  for  the  contour  plot  requests  that  follow.  The 
value  of  n  can  be  from  2  to  20;  default  is  10.  The  values  of  the  xj  can  be 
any  real  number.  The  default  is  a  set  of  contours  equally  spaced  from  the 
minimum  to  the  maximum  in  n  steps.  In  addition,  contour  values  can  be  modi¬ 
fied  by  the  keywords  LOGV  and  LOGD.  LOGV  causes  all  contour  values  to  be  set 
to  values  of  lxlOn  and  3xl0n  from  the  maximum  value  down  for  all  decreasing 
integers  below  n  until  the  number  of  contours  is  exhausted.  LOGD  produces  the 
contours  of  logarithm  base  10  of  the  data.  Contours  are  plotted  by  interpola¬ 
ting  the  data  and  smoothing  to  second  order.  This  can  be  turned  off  by 
including  the  keyword  NOCURV.  The  inclusion  of  the  keyword  COORP  will  cause 
the  zone  size  and  the  values  and  coordinates  of  the  minimum  and  maximum  to  be 
printed  on  each  plot.  If  particles  or  stations  were  included  in  the  calcula¬ 
tion  <N0P>1).  the  particle  positions  can  be  presented  as  dots  on  the  contour 
plots.  This  is  obtained  by  including  the  keyword  PART  in  the  input  deck. 

Histograms  will  be  produced  at  the  mesh  boundary  lower  limit  as  a  default. 
By  specifying: 


PVHST  =  xl,  x2,  dx 

histograms  will  be  produced  along  all  vertical  columns  starting  at  xl  with 
increment  dx  to  x=x2. 

In  order  to  get  the  above  equivalent  plots  from  a  3-D  calculation,  it  is 
necessary  to  specify  XPLANES,  YPLANES,  or  ZPLANES  followed  by  3  values  as 
histograms  were  described  above.  The  following  contour  specifications  will 
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then  be  satisfied  in  those  planes.  If  the  keyword  3DPL0TS  is  present  for  a  3D 
problem,  then  the  user  may  specify  RDCONT,  RECONT,  or  RPCONT  to  plot  a  3D 
representation  01  the  selected  contour  values  when  viewed  in  the  plane  normal 
to  a  line  from  the  mesh  origin  to  the  position  EYE=x,  y,  z.  The  contours  will 
be  considered  opaque  for  this  representation. 

The  prcblem  times  at  which  plots  can  be  selected  are  set  by: 

CTIME=tl,  t2 ,  ...  tn 

DTIME=dt 

FTIME=tl 

LTIME=tn 

If  CTIME  is  present,  it  must  be  the  last  data  in  the  RF.BI.  input  stream.  Plots 
will  be  produced  for  all  times  tl,  t2,  ...»  tn.  The  remaining  three  time 
selectors  will  produce  plots  from  all  dumps  after  tl  at  time  intervals  of  dt 
to  the  last  time  <  tn.  If  dt  is  not  specified,  all  dnmps  between  tl  and  tn 
will  be  plotted.  If  the  keyword  STIME  is  included,  only  standard  times 
(TIMES=1  from  HULL  specification)  will  be  produced. 

B.  JT'TION  PLOTS 


Plots  of  the  station  parameters  can  be  obtained  by  following  the  keyword 
PULL  by  the  keyword  STATION.  Sorting  and  plotting  routines  are  written  to 
file  SAIL.  They  must  be  compiled  separately  and  in  sequence.  The  first 
compilation  must  be  executed  first  to  produce  an  unpacked  station  file  with 
all  time  data  for  station  one  first,  followed  by  all  data  for  station  two 
until  station  NSTN  is  reached.  If  the  HULL  calculation  that  produced  the  file 
was  a  pure  hydrodynamic  run  (STRESS=0),  then  the  execution  of  the  second 
compilation  will  plot  only  the  pressure,  density,  dynamic  pressure,  and 
impulse  as  a  function  of  time.  The  stations  to  be  plotted  can  be  determined 
by  FSTA=nl  LSTA=n2.  All  stations  from  nl  to  n2  will  be  plotted.  The  defaults 
are  nl=l  and  n2=NSTN.  If  STRESS21,  then  a  specific  set  of  plots  can  be 
requested.  These  are: 
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RADIAL  -  produce  plots  of  radial  motion 

AXIAL  -  produce  plots  of  axial  motion 

TOTAL  -  produce  plots  of  total  motion 

STRESS  -  produce  plots  of  each  stress  component 

STRESSDEV  -  produce  plots  of  each  stress  deviator  component 

TENSOR  -  produce  plots  of  the  2nd  stress  invariant  and 
the  principal  stresses 

FLUXM  -  produce  a  plot  of  mass  flux 

We  have  thus  far  not  mentioned  the  use  of  the  problem  number  for  PULL 
input.  It  is  not  needed  unless  it  will  be  used  to  find  the  input  files  from 
the  tape  library  or  the  permanent  file  directory.  If  PROB=XXXXX.YYYY  is 
present,  it  must  be  correct  for  this  run  or  PULL  will  terminate  abnormally. 


SECTION  VII 
RUNNING  A  CALCULATION 


The  preceding  sections  have  addressed  set-up  of  the  various  input  data 
streams  for  running  programs  in  the  HULL  system.  We  have  not  discussed  system 
considerations  because  of  the  frequent  changes  made  in  operating  systems. 
HULL  is  most  often  used  on  CDC  operating  systems.  We  therefore  follow  with  an 
example  control  card  and  data  deck  for  setting  up,  running,  and  plotting  a 
HULL  calculation.  We  will  assume  the  existence  of  HULLIB,  a  user  library 
containing  the  absolutes  of  programs  SAIL  and  PLANK,  and  relocatables  for  the 
HULL  prologue  routines  on  the  HULL  system  file. 


A.  KEEL  TON  DECK 


The  deck  below  will  set  up  50  (X  coordinate)  by  100  (y  coordinate)  zones 
with  DX=DY=.l  cm  from  XQ=0  and  Y0=~5,  respectively.  The  materials  will  be 
air,  iron  (IE),  and  copper  (CU).  Copper  will  be  placed  in  a  cylinder  of 
radius  1  cm  extending  from  Y=0  to  Y=5.  Iron  will  be  placed  in  a  slab  from 
Y=0.5  to  Y=0.  The  iron  slab  is  given  an  initial  velocity  of  2.10*  cm/sec. 

ATTACH  (HULLIB, ID=local) 

LIBRARY  (HULLIB) 

LABEL  (TAPE4, W,T=0,VSN=MINE) 

PLANK. 

ATTACH  (OLD, HULLS YSTEMFILI, ID=local) 

SAIL. 

FTN  ( I=SAIL , B=KEEL , 0PT=2 , L=0 ; 

KEEL. 

*E0R  or  7/8/9  for  card  input  deck 
KEEL  PR0B=1 . 

I MAX =50  JMAX=100  DIMEN=2  STRESS=1 
NM=3  AIR=1  FE=2  CU=3 
REZONEfO  Header 

HULL  EXAMPLE  PROBLEM 
MESH  X0=0  XMAX=5  Y0=-5  YMAX=5 
GENERATE 

PACKAGE  AIR  RECTANGLE 

DELETE  RECTANGLE  XR=1  YB=0 
RECTANGLE  YB=0.5  YB=0 
PACKAGE  CU  RECTANGLE  XR=1  YB=0 
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PACKAGE  FE  V=2.E5 

RECTANGLE  YB=-5  YT=0 

END 

•EOR 

SAIL  LINENO 

Change  cards  needed  by  this  run 

B.  HULL  1DN  DECK 


ATTACH  (HULLIB,  II>=  local) 

LIBRARY  (HULL IB) 

LABEL  (TAPE4 , R, VSN=MINE) 

PLANK. 

ATTACH  (OLD.HULLSYSTEMFILE, ID=local) 
SAIL. 

FIN  ( 1= SAIL, B= HULL, OPT=2 , L=0 ) 

HULL. 

*EOR 

SAIL  LINENO 
local  changes 
•EOR 

HULL  PR0B=1 .  CYCLE=0 
INPUT 

TIMES=3  DMPINT=5 . E-6 

PTSTOP=50.E-6 

MATPROP 

MAT=2  YLDMAX=15.E9 
MAT=3  YLDMAX=4.5E9 
ENDPROP 
•EOR 


C.  POLL  SUN  DECK 


ATTACH  (HULLIB, ID=local) 

LIBRARY  (HULLIB) 

LABEL  (TAPE4,R.VSN=MINE) 

PLANK. 

ATTACH  (OLD.HULLSYSTEMFILE, ID=local) 

SAIL. 

FTN  ( I=SAIL, B=PULL, L=0 , OPT=2 ) 

PULL. 

•EOR 

DCONT=7 . 8  7.9  8  8.1  8.2  8.3  8.6  8.7  8.8  9 

FVHST  WE  CT  COORD  NOCURV 

STIME 
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SECTION  VIII 
SENSE  SWITCH  CONTROL 


HULL  can  be  controlled  for  set  up  by  SAIL  and  during  running  by  CDC  sense 
switch  settings.  If  the  sense  switches  are  set  before  PLANK  runs,  they  have 
the  following  effects: 

SENSE  SWITCH  ON  BEFORE  PI.ANK  RUNS 

SENSE  SWITCH  NUMBER  EFFECT 

0  (All  switches  off)  ECS  store  of  the  mesh 

5  In-core  central  memory 

storage  of  the  mesh 

6  Disk  storage  of  the  mesh 


If  the  sense  switches  are  turned  on  after  HULL  or  PULL  have  started,  the 
following  actions  will  be  taken: 

SENSE  SWITCH  ON  AFTER  HULL  OR  PULL  START 

SENSE  SWITCH  NUMBER  EFFECT 

1  Stop  at  end  of  this  cycle  (HULL) 
plot  this  dump  and  stop  (PULL) 

2  Route  current  output  to  printer  (both) 

3  Display  current  time,  cycle,  tire  step 
and  estimate  of  running  time  on 
operators  console  (HULL) 

Display  time,  cycle,  and  type  of  plot 
being  made  (PULL) 

4  Change  dump  tapes  (HULL) 

Change  to  next  input  tape  (PULL) 

5  Start  double  buffering  of  a  disk  storage 
problem.  Increases  central  memory  but 
develops  I/O.  'HULL) 
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5  OFF  Stop  double  buffering,  return  unneeded 

central  memory.  (HULL) 

6  Roll  ECS  out,  resume  by  operator  GO.  (HU1L) 

The  sense  switches  can  be  toggled  on  by  the  operator  at  his  console  or 
through  the  use  of  control  cards  in  the  input  stream. 
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Ammon  a 

BASIC  BCD AXIOMS  AND  SOLUTION  ALBOftTIHM 


The  HDLL  code  solves  the  partial  differential  equations  of  continuum 
mechanics.  Explicit  terms  for  heat  conduction  and  viscous  effects  are  not 
included.  The  equations  solved  in  finite  difference  form  by  HDLL  are: 


p  +  pu.  .  =0 

1.1 


(A-l) 


P*j  -  Tij.i  =  “ 


(A-2) 


where: 


PE  (Tijuj)'i  = 

Tij  =  f(P'X'ui.j  +  Uj,i> 


(A— 3 ) 
(A-4) 


=  material  density  (gm/cc) 

=  material  velocity  (cm/sec) 

=  SAj  -  stress  tensor, 

P  =  hydrostatic  pressure  (dynes/cm*) 

=  stress  deviutor  (dynes/cm*) 

=  gravitational  body  force  (cm/ sec*) 

=  I  +  1/2  ujUj ,  total  specific  energy, 
I  =  internal  specific  energy  (ergs/gm) 


subject  to:  S ^ j  S A .  X  =  flow  stress. 


The  devia'.^ric  strain  rate  tensor  can  be  expressed  as: 


'ij  =  *ij  “  <ui,j  +  °j,i>/2  -  6ij"i,i/3’ 


(A-5) 


The  deviatoric  stress  tensor  is  linearly  related  to  the  strain  tensor  by: 


Si.  =  20.^ 

where  G  is  the  elastic  shear  modulus. 


Plasticity  is  modeled  by  the  Von  Mises  yield  criterion: 


S, .  S-.  <  Y  =  tensile  flow  stress. 

ij  3 


(A-6) 


If  the  stress  deviators  are  not  contained  within  this  measure,  they  are 
brought  back  to  the  yield  surface  by: 


S. . 

ij 


=  S. 


ij 


JlL 


,1/2 


3SijSij 


(A-7) 


The  initiation  of  failure  can  be  determined  by  one  of  three  criteria: 

(1)  Maximum  principal  stress  can  be  used  as  a  crude  means  to  deter¬ 
mine  if  a  spall  threshold  or  failure  in  plane  strain  has  occurred.  The 
principal  stresses  are  defined  by  solving: 

|  =  0.  (A-8) 

Failure  is  initiated  whenever: 

Maximrm  [a]2.  F’Y 

where  F  is  a  constant  between  2.5  and  4. 

(2)  Maximum  principal  strain  can  also  be  used  to  determine  the 
ouset  of  failure  in  plane  stress  or  ductile  failure  by  solving  the  analog  of 
Equation  (A-8)  for  the  maximum  principal  strain  as: 


and  assuming  failure  occurs  whenever 
Maximum  [X]  2  ef 

where  e^  is  the  tensile  strain  to  failure  of  the  specific  material. 


(A-9) 
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(3)  Triaxial  states  at  which  failure  can  occur  are  based  on  a  P/Y 
model  (Reference  6).  In  this  formulation  failure  is  initiated  when  the  maxi¬ 
mum  principal  strain  (X,  Equation  (A-9))  exceeds  a  failure  surface  =  f(P/Y) 
which  is  determined  by  experiment.  Failure  in  HULL  is  initiated  by  inserting 
a  minimum  numerically  significant  void  in  the  failed  zone  and  allowing  the 
tensile  forces  to  relax.  Subsequent  tension  or  shear  of  this  region  is  not 
supported  (i.e.,  the  zone  defines  a  slip  surface  and  tensile  stresses  are 
zeroed) . 


Rewriting  Equations  (A-l)  through  (A-3)  in  axisymmetric  cylindrical  coordi¬ 
nates; 


d(ru) 
~Ji  + 


3v 

31 


]-• 


(A-10) 


pft  - 


(A-ll) 


a  T  dT 

p*  "  -JFT  “  ~ST  =  ~pg  (A-12) 

PE  '  5T  [  t(»Trr  +  'Trz>  ] 


Tz  UTrz  +  vTzz>  =  "Pv8 

r,z  =  radial  and  axial  coordinates 
u  =  radial  velocity  component 
v  -  axial  velocity  component 


T 

Arr 

= 

Srr-P 

T 

zz 

= 

sz2  -  P 

Tee 

= 

s«e  “  p 

T 

rz 

2: 

Tzr  =  "S 

(A- 13) 
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rr  *  '  D) 


S‘„  -  20(J„  -  D) 

SSe  ■  -  ^ir  "  &ZZ 
St.  -  20<L.) 


.  du  dv  . 

rr  ~  <Tr  *  8zz  “  Tz 


_  r  dru  .  Jv  v 

D  =  l  5J*  +  3TS  j/3 


Subject  to  the  yield  criteria: 


SaH  "  S^'  J  ’  «}r  *  »},  -  S>,  *  S„S„)  <  3 


X* 


a  1/2  * 

s«p  [  jy]  J  >  ]r 


Failure  implies  that  *£«frac,  where: 


T„-T, 


=  1/2(T„  +  Tzz>  +  [  (-^V^)  +  Trz2  ] 


1/2 


wfrac  *s  1*0*11*  yield  stress,  a  number  usually  between  2  and  5  times 
the  flow  stress. 

D1FFKSEHCR  4ETMD  (20) 

Equations  (A-1Q),  (A-ll),  (A-12),  and  (A-13)  are  solved  on  a  finite 
difference  mesh  of  discrete  intervals  Ar^  ,Ai-  in  the  radial  and  axial 
coordinates,  respectively.  The  solution  is  advanced  explicitly  from  the 
initial  conditions  by  discrete  time  steps  Atn.  The  solution  is  thus  defined 
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on  the  mesh  (r^,  Zj 
space  defined  by 


tn)  where  each  of  the  variables  f(r,z,t)  in  the  solution 


fCrj,  Zj,  tn) 


where 


r . 
1 


i-1  Ar . 

*.  *  ll  **  J ♦ -r 


"k=l 


(A-14) 


z  . 
J 


2 


+ 


Az . 
_ J 

2 


(A-15) 


n 

t“  -  tc  +  j  Atn  (A-16) 

k=l 

It  should  be  noted  that  this  subscript  notation  is  distinct  from  the  previous 
tensor  notation. 


V  2o'  *o  are  the  points  in  the  r,  z,  and  t  coordinates,  respec¬ 

tively.  We  also  find  it  convenient  to  define  Ari+1/2  =  ti+1  -  xi  and  Azj+1/2 
Zj+1  "  zj* 

The  variables  are  defined  for  each  (r.,  z^)  at  time  tE*  Spatial 

derivatives  will  be  approximated  by  the  finite  difference  operators  6r  and  8z 

as: 


6r 

6Z 


(fi.j)  "  <*?♦!/' .J 


(fn  )  =  (fn 

Ui.j+l/2 


-  fi-m.jJMr, 
"  fi,j-l/2)/Azj 


5 


Time  derivatives  will  be  approximated  by  (f ^  j  -  f?j)/At  where  At  =  Atn  is 
assumed.  It  should  be  noted  that  the  time  derivatives  appearing  in  Equations 
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(A-10)  through  (A-16)  are  Lagrangian;  that  is,  they  refer  to  a  coordinate 
system  moving  with  the  material.  Hence  the  fi  ,  values  refer  to  the 
property  f  at  some  displaced  location.  This  must  be  borne  in  mind,  as  it  is 
not  noted  explicitly  in  the  equations  which  follow: 

The  solution  for  (A-ll),  (A-12),  and  (A-13)  is  thus: 

»”+1  »„■>-*£[  -  l6'rs?r-6«s»2  +  Sjdk  ]  (4-17) 


sZ(sL  -  P“"1/2)  ] 


(A-18) 


En+1  =  En  +  At  f  1  rjrrun+l/2^sn  _  pn+l/2)  +  rvn+l/2sn  , 
„n  L  i  rr  r  rzJ 


(A-19) 


6z[vn+l/2(sn  _  pn}  +  nn+l/2gn  -  1 
zz  rxJ  J 


Since  the  pressure  is  a  function  of  density  and  internal  energy  only,  we 
can  obtain  a  substantive  time  derivative  for  pressure  by 


1  3P  .  3P  • 
P  ~  ^  P  31  1 


(A-20) 


From  Equations  (A~ll)  through  (A-13)  we  can  obtain: 


Pl  +  ^1.1=  0 


(A-21) 


Combining  Equations  (A-l) ,  iA-20) ,  and  (A-21)  we  can  write 


;  _  r  ap  ^  p  9P  1 
p  t  3p  p2  31  ^ 


bat  the  term  in  parentheses  is  the  square  of  the  isentropic  sound  speed 
hence : 


P  = 


i»  i 


(A-22) 
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Using  Equation  (A-22)  we  can  calculate 


pn+1/2 

1+1/2, j 


_  pu  _  (Dr  a)n 

"  *i+l/2 ,j  2  s  'i+1/2 ,j 


(ru) 


i+1/2,  j 


(A-23a) 


00+1/2 

i, j+1/2 


"  Pi, j+1/2 


At 

T 


(pcs*) 


i, j+1/2 


&Z  (vi,j+l/2> 


(A-23b) 


where : 


pn 

i+1/2, j 


Pi+l,j  P?,j  +  p?,j  >M,j 

pi,j  +  p°+l.j 


(A-24a) 


pn 

i, j+1/2 


pn  n 

^.j+l  p 


+  pn  or 

i,j  Pi«j  pi.j+l 


+ 


P 


n 

i.  j+1 


(A-24b) 


(pCS'>;+l/2,j  -  Mi»  [  <pO?,j  !  ] 

(oC.*,!.3*1/2  ■  Mi”  [  «’CS‘)i.j  i  ] 

2  -j 

and  (pC  ) -  •  is  obtained  from  the  equation-of-state, 
S  1#  J 


(A-25a) 

(A-25b) 


The  spatial  centered  pressures  of  Equation  (A-24)  were  chosen  to  avoid 

over-acceleration  in  regions  of  large  density  gradients.  This  formulation 

results  in  correct  free  surface  conditions  at  solid/air  interfaces  at  one 

limit  and  degenerates  to  a  simple  average  at  the  other  limit  of  small  density 

2 

variations.  Minimum  values  of  pCs  in  Equation  (A-25)  are  used  to  keep  from 
overpredicting  the  pressures  calculated  in  Equation  (A-23).  The  physical 
rationale  for  this  choice  is  the  assumption  that  the  most  easily  compressed 
material  will  restrict  the  pressure  rise  since  it  will  compress  first.  The 
procedures  also  switch  the  scheme  from  second  order  accuracy,  in  regions  of 
smooth  flow,  to  first  order  accuracy  in  regions  of  rapid  transitions  in  a  more 
simplified  calculational  methodology  than  described  by  Harten  and  Zwas 
(Reference  3). 
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Time  advanced  velocities  for  Equation  (A-19)  are  obtained  by  differencing 
Equations  (A-14)  and  (A-15)  by: 


„n+l/2  _  n  At  K,/pn  v 

i+1/2, j  "  ui+l/2, j  ,*+1/2  '  1+1/2, j' 

Pi+l/2,j 


(A-26) 


n+1/2  _  n  _  <■  ,pn  . 

i, j+1/2  vi, j+1/2  ~  ,  n+1/2  i, j+1/2' 

Pi, i+1/2 


At 

2  ® j+1/2 


where : 


8j  +  «j+l 


8 j+1/2 


i+1/2. j 


nn  n  .  n  n 

Pi» j  gl» j  +  pi+l.j  pi+l,j 

Pi,j  +  pi+l,j 


(A-27) 


i, j+1/2 


„n  n  n  n 

Pi,j  vi,j  *  pi, j+1  vi, j+1 
Pi. j  +  Pi+l,j 


Pi+l/2, j  "  ***  <pi.j  '  pi+l,j> 


(A-28) 


Pi, j+1/2  "  <pi,j  '  pi,j+l) 

where  again  the  selection  of  spatially  centered  quantities  in  Equations  (A-27) 
and  (A-28)  was  made  to  reduce  overacceleration  in  regions  of  large  density 
gradients. 


The  density  Equation  (A-10)  is  also  solved  in  the  first  step.  The  new 
densities  are  stored  and  used  in  the  fluxing  (or  Rezone)  calculation  below. 
Equation  (A-10)  is  differenced  as: 


p*+1.  =  p?  .  \  l  -  —  6r  <ru)?+V2  -  At  bzy^1/2  ] 
i,j  pi,j  l  A  r.  i.J  i,j  J 


(A-29) 
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using  the  definitions  above  for  the  time-centered  velocities. 

The  deviators  by  definition  have  the  property  =  0  leaving  only  Srr> 
S  and  S__  to  be  retained  as  stored  quantities  from  the  previous  time  step. 

I*  ZZ 

These  deviators  are  defined  at  each  mesh  point  in  the  solution  space.  The 
numerical  solution  is  obtained  explicitly  by: 


»  P°.W  *  >*.j> 

S“>i,jtl/2  -  ®i,j+l/2  f  ♦  Pi ,j+1 


»i.l/2,j  “  <VFi.j  •  VFi+i<J) 

®i,j+l/2  ”  Min  *FFi,j  '  FFi,j+l* 

VF •  .  =  Fractional  volume  of  solid  in  zone  i,j. 

J 

Other  quantities  required  above  are  defined  by  Equations  (A-17),  (A-18),  (A- 
19).  and  (A-27). 

At  the  end  of  this  step,  the  stress  deviators  themselves  are  updated  by: 


=  Srr“  .  *  20  it  [*'  “  »] 

a  1  a  a  •* 


rr 


(A-34a) 


i.J  i.J 


Szz”+1  -  S„°  .  +  20  At  [0Z  -  D] 

i.J  i.J 


(A-34b) 


Srz”+1  =  Srzn  +  G  At  (8r  (v)i  i  +  8z  <n>i  i> 

“i.J  “i.j  1#J 


(A-34c) 
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Where  G  is  a  material  dependent  modulus  and: 


O  .  r  il  (rn)fV2  +  6,  v?+1/2  1/3 

*•  rj  i.j  i.J  J 

These  definitions  of  the  stress  deviators  are  subject  to  the  yield  criteria:. 

S  a  =  Sea  ;  J  =  (S2  +  S2  +  S2  +  S  S  )  <  1/3  Y2 
a0  a0  rr  rz  zz  rr  zz  ~ 

2  1/2 

*  Sa0  <«->  J  >  1/3  l2 


Failure  implies  that  °A®frac  ,  .rhere: 


o  =  1/2 


(Trr  + 


Tzz> 


T  -  T 


+  T 


rz 


and  ®£rac  is  the  tensile  yield  stress,  a  number  usually  between  2  and  5  used 
in  the  HULL  code.  The  modularity  of  the  code  makes  substitution  by  more 
comprehensive  failure  models  an  easy  matter. 


Implementation  of  failure  in  a  zone  is  accomplished  by  inserting  a  void 
or  gas  (air)  into  the  “unoccupied*  volume  of  a  cell  after  allowing  the  failed 
solid  to  relax  to  ambient  density.  Use  of  a  gas  rether  than  a  void  provides  a 
damping  mechanism  reducing  the  tendency  of  a  failed  cell  to  “chatter*  as  it  is 
reloaded. 

The  definition  of  a  slip  surface  is  the  existence  of  more  than  one  solid 
(but  no  gas)  in  a  cell.  Such  a  cell  will  permit  compressive  stress  but  will 
not  support  shear  or  tensile  force.  Cells  which  contain  a  gas  or  void  will 
not  support  either  distortional  or  spherical  tensile  forces.  This  implies 
that  surfaces  across  failed  regions  are  automatically  defined  as  slip  surfaces 
whenever  a  failure  criterion  is  invoked. 
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•itAr 


D1FFEUBNCZ  MBHOD  (3D) 


Equations  (A~l)  through  (A-3)  are  solved  on  a  finite  difference  mesh 
composed  of  discrete  spatial  intervals  xj,  y j,  z^  the  three-dimensional 
Cartesian  coordinate  system  having  coordinates  x,  y,  and  z.  The  solution  is 
advanced  explicitly  from  the  initial  conditions  by  discrete  time  intervals 
Atn.  The  solution  at  any  time  tn  for  any  variable  f(x,y,z,t)  in  the  solution 

space  is  n  defined  by  f^k  =  f(xj,yj,z^,tn)  where: 

i-1 

\  Ax. 

Xi  =  xo  +  1  +  T1 


yj  =  y0  + 


*  % 


'k  ■  zo  ♦  1  ?*. 


♦T* 


*”  ■  lo  +  l  A*. 

m=l 

Note  that  on  this  mesh  Axj  =  xj+i  -  x^» 

This  definition  holds  analogously  for  y,  z,  and  t. 

The  variables  that  are  defined  for  each  point  in  the  mesh  at  time  tn  and 
location  x.,  yjj  Zk  are; 

P  =  pressure  {dynes/cm* ) 
u  =  x  component  of  velocity  (ca/sec) 

v  =  y  component  of  velocity  (cm/sec) 

w  =  z  component  of  velocity  (cra/sec) 

I  =  internal  specific  energy  (ergs/gm) 

M  =  mass  (gms) 

Mn  =  mass  of  the  nth  species  (gms) 

\*  =  volume  occupied  by  the  nth  species  (cm*) 
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In  =  total  internal  energy  of  the  n**1  species  (ergs) 

Sxx  =  a  coordinate  deviatoric  stress  (dynes/cm1) 

Syy  =  y  coordinate  deviatoric  stress  (dynes/cm*) 

S  =  shear  stress  nomal  to  z  coordinate  (dynes/cm*) 
xy 

sxz  =  shear  stress  normal  to  y  coordinate  (dynes/cm*) 

S  =  shear  stress  normal  to  x  coordinate  (dynes/cma) 
vz 

Since  the  deviators  have  the  property 
,<?zz  _  _Sxx  _Syy 

S  is  not  explicitly  stored  as  a  mesh  variable, 
z  z 

If  ire  define  the  central  difference  operator  6x(f.  .  ,  )  as: 

8X(£i,j,k)  =  ( fitiy.2,1iik_~.  ,f  i-i/ji^ 

^i 

then  we  can  express  the  calculation  for  the  velocity  components  as: 

nn+l  =  nn  -»  [  6x(sn  _  pn+l/2)  +  8y(Sn  }  6z(Sn  }  1 

L  xx  xy'  °  xz'  J 

V”+1  =  +  ^  [  6y(S^y  '  P”+1/2>  +  SX(S^y>  +  &Z(SyZ>  ] 

wn+l  =  w»  +  ^  [  6Z(-S*x  -  -  P“+1/2>  +  &X<SZZ)  +  6y(Syz)  ] 

where : 

*1+1/2  =  Pi+l/2  ~  T  pCi+l/2  6X(ni+l/2) 

and: 


n  Pj+lPi  PiPi+l 

i+1/2  Pi  +  Pi+1 
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<1/2  -  MN 


K.j.h 


tti+l/2 


piui  +  Pi+ltti+l 


+  pi+l 

Solution  o£  the  conservation  of  energy  equation  is: 


En+1  =  En  +  At  r  x  (_p  un+l/2  +  sn  Qn+l/2  +  sn  yn  +  §n  „n) 
n  L  xx  xy  xz 

P 


+  &y  (-pyU+1/2  +  sn  vn+l/2  +  sn  un  +  sn  wn) 

yy  xy  yz 


6Z  (-Pwn+1/2  4  gn  wn+l/2  +  Sn  u®  +  Sn  v®) 
zz  xz  yz 


where  the  boundary  value  velocities  at  time  n+1/2  are  given  by: 


n+1/2  n  At  tr,  n 

ui+l/2  “  ui+l/2  7  n+1/2  (Pi+l/2 

2pi+l/2 


with  similar  forms  for  v  and  w  with 


Vl/2  = 


PjUj  P  j+1  ®j+l 
Pi  +  Pi+i 


n+l/2_  n  -  At  .x  n  . 

piH/2  ”  pi+l/2  '  "  2  ui+l/2; 


At 


n  n  n 

Pi+l/2  =  Maximum  (pi  ,  pi+1). 


Each  of  the  SQp  terms  are  calculated  in  a  manner  similar  to  SM  as 

-W  X  *  Pi  s„.+1 


Pi  +  Pi+1 


*i+l/2  =  MINtvfj,  vfi+1) 
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and  vf^  =  the  fraction  of  zone  i,  j,  k  which  is  occupied  by  a  solid  material. 
The  time  n+1  deviators  are  calculated  by: 

Sxx1  =  Sxx  +  2G(6x(nn+1|/2)  -  D)At 

S^1  =  +  2G(8y(vn+1/2)  -  D)At 

S”^1  =  S^.  +  G(6x(vn)  +  6y(un)>At 

s”*1  =  szz  +  G(8x(»E)  +  6z(un;)At 

Sy^1  -  SyZ  +  G(8y(w11)  +  ftz(vn))At 

where : 

D  =  I6*(un+1/2)  +  Sy^1/2)  +  8z(w®+1/2)]/3. 


The  stresses  are  then  compared  with  the  yield  surface  definition  by: 


J  =  S2  +  S2  +  S  s  +  sa  +  S2  +  S2  <  Y2/3 

xx  yy  xx  °yy  °xy  axz  ’’yz  -  x 


If  this  criterion  is  not  met,  then  each  stress  component  is  reset  by: 
r  Y2  -il/2 

Sa0  L  3 J  J  SaP 

where  represents  each  of  the  elements  of  the  stress  tensor. 

ADVKCTION  PHASE  CALCULATION 

Up  to  this  point  the  entire  calculation  has  been  in  a  Lagrangirn  frame¬ 
work.  The  next  step  is  to  indicate  the  volume  change  due  to  spherical  forces 
or  pressure  and  prepare  the  system  for  the  fluxing  calculation.  This  is  done 
by  solving  the  continuity  equation  by: 

vn+! 


=  Vn(l+3AtD) 


where  V  is  the  cell  volume.  The  equation  of  state  is  then  called  and  the 
pressure  in  the  new  zone  is  calculated.  If  the  zone  contains  more  than  one 
material,  an  iteration  procedure  is  used  to  equilibrate  the  pressure  of  each 
material  species  in  the  zone  to  a  common  value.  This  is  dons  by  reparti¬ 
tioning  the  zone  among  the  materials  present.  Species  energy  is  calculated  by 
applying  At  =  P  V  as  the  volume  V  occupied  by  each  species  at  pressure  P 
subject  to  conserving  the  total  mass,  volume,  and  internal  energy  of  the  zone. 


The  calculation  is  returned  to  a  fixed  (or  uniformly  translating  or 
expanding)  grid  by  reapportioning  material  from  the  displaced  to  the  fixed 
grid.  The  volume  to  be  fluxed  is  defined  as  the  sum  of  the  transfer  volumes 
in  each  of  the  three  directions,  which  are  given  by: 

Vi+l/2  =  (xi+l/2  "  Xi+1/2J  Ayj  Azk  ^  1/V° 


Vj+l/2  ~  (yj+l/2  A+lll)  Asi  Azk  V°  1/VQ 
Vk+l/2  =  (zk+l/2  “  Zk+1 / 2 '  Axs  Ayj  ^  1/Vn 


where 


Vi+l/2  =  voluae  (ia  cm*)  to  be  transported 

between  cell  i,j,k  and  cell  i+l,j,k 


(A-35) 


(A-36) 

(A-37) 


n+1  c  vn  .  „n+l  *«■ 

‘i+1/2  xi+l/2  it?/?  At 


n 

x-  + 


i+1/2 


xi-: 


n+1  n  n+1  n 

Uj  pi  +  ui+1  pi+1 

pi  +  pi+l 

with  analogous  definitions  for  j  and  k. 

The  transfer  volume  defined  by  Equations  (A-35)  -  (A-37)  is  apportioned 
among  the  available  materials  in  the  donor  cell  in  accordance  with  the 
diffusion  limiter  algorithm.  This  algorithm  is  ad  hoc  but  has  been  formulated 


n+1 

i+1/2 


u 
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over  the  years  by  observation  of  various  numerical  experiments  and  has  been 
found  to  be  modestly  successful.  The  algorithm  uses  nearest  neighbor  informa¬ 
tion  in  an  attempt  to  unmix  mixed  material  cells.  This  is  done  by  assigning  a 
relative  transport  weight  to  each  material  in  the  donor  cell,  and  then 
normalizing  to  the  transfer  volume  given  in  Equations  (A-35)  -  (A-37).  Each 
direction  (up  to  6  in  3D)  is  handled  separately  subject  to  the  constraint  that 
the  sum  of  the  exiting  volumes  for  each  material  shall  not  exceed  the  initial 
volume  of  that  material. 


The  relative  transport  weight  W  for  each  material  n  is  defined  by 
VnK  -  VnU 


w  =  r  ^  + 1.0 1 

n  1  VnE  +  Vn0  J 


=  [ 


2V 


nR 


(A-38) 


VnR  +  VnD  J 

where  =  fractional  volume  of  material  n  in  the  receiver  cell 
=  fractional  volume  of  material  n  in  the  upstream  cell. 

The  fractional  volume  is  defined  as  the  volume  of  material  n  divided  by  the 
total  cell  volume.  Two  exceptions  to  Equation  (A-38)  are  noted.  If  there  is 
no  material  n  in  the  donor  cell  then  Wn  =  o,  and  if  VnR  and  VnD  are  zero  then 
V  =1. 

Finally  the  transport  fraction,  for  each  material  n  in  a  particular 

direction  m  is  obtained  by  normalizing  the  relative  transport  weights  by 


pm 


pnr* » 

L  n  n 


where  Vp  =  transfer  volume  given  by  Equation  (A-35),(A-36j,  or  (A-37) 

v  ~  volume  of  material  n  in  donor  cell. 

Tn 

am 


The  transport  fractions  A 


are  subject  to  the  constraint  that  the  sum  over 


all  directions  (six  in  3D)  shall  not  exceed  unity  for  each  material  n.  If  it 
does,  snitable  corrections  are  mado. 

The  transport  fraction  applies  not  only  to  the  volume  but  also  to  the 
mass  and  internal  energy  for  material  n.  That  is,  the  fraction  Anffl  of  each  of 
these  quantities  will  be  transferred  from  the  donor  to  the  receivor  cell.  The 
sum  over  n  for  each  direction  defines  the  total  flux  of  mass  and  internal 
energy. 

The  transport  methodology  will  be  illustrated  for  one  direction  (flux  to 
the  right).  The  remaining  directions  are  analogous.  Consider  cell  i»j,k  at 
the  end  of  the  Lagrangian  phase,  denoted  above  as  time  n+1.  For  each  material 
n  the  amount  of  mass  transporting  is 

6*£  =  (AnR  «n) donor 
while  the  volume  transporting  is 

&VR  =  (A  „  v  )j 
n  nR  Yn' donor 

and  the  internal  energy  is 

6(MI)n  =  tAuR<MI>n]donor- 

The  total  mass  fluxing  (to  the  right)  is  therefore 

6#  =  \  6# 

L  n 

c 

while  the  internal  energy  is 

6(MI)R  =  £  6(Mnin)R. 


The  total  energy  is 

6(ME)R  =  6(MI)r  +  1/2  SM®  (Uji  +  v<Ja  +  Wd») 
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where  d  indicates  donor  cell  velocities.  The  transporting  material  carries 
with  it  a  proportionate  amonnt  of  momentum,  thus 

6(Mu)®  =  n^6M® 

6{Mv)S  =  v.b*^ 
a 

6(Mw)R  =  w^&iJ8 

With  most  options  HULL  transports  certain  histories  from  cell  to  cell.  These 
histories  essentially  transport  with  the  mass;  hence,  for  the  kth  history  H 

»  Hdt 

where  d  again  indicates  the  donor  cell. 

Similar  terms  are  obtained  for  the  other  directions  fore  (F)  and  above 
(A),  and  from  previous  calculations  for  left  (L),  aft  (T),  and  below  (B). 

Terms  such  as  the  mass  M.n+1  v  are  redefined  from  the  displaced 

Lagrangian  mesh  to  the  original  fixed  (or  uniformly  translating)  Eulerian  mesh 

eM.n+l  as  follows.  For  material  n 
i»J»k 

eMn+1  =  HP+i  +  6#  +  6Mt  +  6#  -  6Kr  -  6#  -  6mA 
r  n  nnnnnn 

eV»+l  =  Vn+1  +  6VL  +  6VT  +  6V®  -  6VR  -  bV8  -  bVA 
n  n  nnnnnn 

e.'MI)a+1  =  (MI)n+i  +  6(MI)L  +  6(MI)t  +  6(MI)b  -  b(MI)8  -  b(MI)F  -  b(MI)A 
n  n  nnnnnn 

The  velocity  is  defined  by  conserving  momentum  with 

eua+1  =  [un+1  Mn+1  +  6 (Mu) L  +  6 (Mu)T  +  6(Mu)B 

-  6(Mu)R  -  6(Mu)F  -  6(Mu)A]  /£  e*P+1 

n 


evn' 


Hi  _ 


[vn+1  M1  *  +  b(Mv)L  +  6 (Mv)T  +  6(Mv)B  -  b(Mv)R 


-6(Mv)f  -  6(Mv'Al  /  J  ®M“+1 

n 

cwn+1  =  [wn+1  Mn+1  +  6(Mw)l  +  6(Mw)T  +  6(Mw)B  -  6(Mw)S 

-6(Mw)F  -  6(Mw)A]  /  ^  °Mn+1 

L.  n 


The  new  internal  energy  is 

CIn+1  =  [In+1  Mn+1  +  6(ME)L  +  6(ME)T  +  6(ME)B  -  6(ME)R  -  6(ME)F  -  6(ME)A] 

/  \  eMn+1  -  1/2  [(eun+1)2  +  (evn+1)*  +  (ewn+1)*] 
l  n 

n 


The  new  histories  are  defined  by 

eH“+1  -  (H£+1  Mn+1  +  6h£  +  6H?  +  6H|  -  6H®  -  61^  -  6H£)/  £  ejpi+1 


It  should  be  noted  that,  in  general, 

(eMn+1)(eIn+1)  it  \  e(MI)®+1 
n 

as  these  terms  are  defined  above.  This  is  due  to  a  dissipation  of  kinetic 
into  internal  energy.  Total  energy  is  conserved  by  adjusting  the  individual 
energies. 


APPENDIX  B 
HOUATIONH-OF-STATE 


The  equations-of-state  do  not  generally  attempt  to  aconrately  cover  all 
regions  of  P,  V,  I  space.  They  were  developed  to  provide  valid  material 
response  in  problems  of  interest  in  conventional  munitions  research  -  shaped 
charge  jet  formation,  self  forging  fragment  formation,  metal-on-raetal  penetra¬ 
tions,  bomb  penetrations  into  concrete,  and  other  similar  problems.  This 
emphasis  has  meant,  for  example,  that  great  care  was  given  to  modelling  solid 
and  liquid  metal  behavior,  but  metal  vapor  behavior  was  simpl  fied  to  an  ideal 
gas  formnlation.  Problems  of  conventional  munitions  research  interest  do  not 
usually  require  highly  accurate  metal  vapor  response. 

The  HULL  equations-of-state  (EOS  6)  have  been  implemented  in  EP1C3  with 
some  minor  simplifications.  Thus,  this  report  can  serve  as  documentation  for 
many  of  the  EPIC3  equations-of-state  as  well  as  those  in  HULL. 

In  what  follows,  the  term  '“equations-of-state.11  will  refer  to  the  determina¬ 
tion  of  all  stress-strain  relations  including  deviatoric  as  well  as  hydrosta¬ 
tic  behavior. 

1.  Metal  Eqaation*-of-State 

The  metal  equation-of- state  relationships  are  based  on  Mie-Gruneisen  and 
gamma-law  gas  pressure  formulations  and  the  von  Mises  yield  criterion.  The 
yield  strength  can  be  a  bi-linear  function  of  accumulated  plastic  strain  (work 
hardening)  and  a  tri-linear  function  of  internal  energy  (thermal  softening). 
The  internal  energy  density  (energy  per  unit  mass)  at  onset  of  melt  and 
vaporization  are  increased  with  increasing  density.  Increasing  internal 
energy  is  not  allowed  to  increase  pressure  during  the  melting  and  vaporization 
ph  ase  changes.  We  have  made  no  provision  for  solid-solid  phase  changes  al¬ 
though  such  phase  changes  could  be  added  relatively  simply. 
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A.  The  Pressure  Equations 


The  basic  HULL  pressure  equation  for  a  solid  or  liquid  metal  is  the 
Mie-Grunoisen  equation: 

P  =  PH  <l-ni/2)  +  p(I-I0)  +  P0 

where : 


Pg  is  the  Hugoniot  pressure  given  by: 
PH  =  Bxp  +  Ba;2  +  B,;»  if  JJ>0 

=  B^p  if  P<0 


and  p  is  the  excess  compression: 

P  =  p/p0  "  1 

p  is  current  density 
p  is  ambient  density 

O 

dP 

f is  the  Gruneisen  parameter  (V  8iven  by: 


r  -  rop0/p 

I  is  the  internal  energy  per  unit  mass 

10  is  the  internal  energy  per  unit  mass  at  ambient 
density  and  pressure  (1  atmosphere) 

Ps  is  ambient  pressure  and  is  calculated  from 

,1.013  x  10*  dynes/cm2.  ...  ,,  „  * 

P0  =  ( - j - 1 - ).  Min  (I,  I0) 

Ao 

This  allows  Pe  to  be  1  atmosphere  when  I  =  I0  and  to  decrease  smoothly  to  zero 
as  I  drops  below  the  ambient  level  (I0)» 


The  Hugoniot  coefficients  are  calculated  from: 


Bx  *  P0Co* 

Ba  =  Bx  [1+2  (s-1)] 

B,  =  Ba  [2  (s-1)  +  3  (s-1)*] 


where  Cq  and  s  define  the  linear  shock  velocity  (Us)  versus  particle  velocity 

(U  )  relationship. 

P 


The  assumption  of  a  linear  U  versus  U 

s  p 


relationship: 


U  =  C  +  s 
s  o 


combined  with  the  Hugoniot  jump  conditions  actually  results  in  the  Hugoniot 
pressure  equation: 


PH  =  poCo*  +  1)/U  '  W(s-l)]1 


Expanding  the  denominator  by  the  binomial  theorem,  simplifying  and  deleting 
all  terms  with  the  multipliers  beyond  p*  yields  the  coeffients  Ba,  Ba  and  B». 
This  expansion  process  has  also  deleted  the  pole  in  the  formulation. 


The  Gruneisen  parameter  assumption,  namely 


r  ■ 

is  considered  reasonable  for  metals  (Reference  1)  and  results  in  well-behaved 
pressures  at  high  values  of  density.  The  other  often  used  assumption  that 
does  not  change  with  increasing  density  results  in  pressure  maxima  for  rela¬ 
tively  small  density  increases. 


The  internal  energy  density,  I,  used  in  the  Nie-Gruneisen  pressure 
calculation  is  not  always  the  zone's  internal  energy  density,  1^.  The  rela¬ 
tionship  is  seen  below: 
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where  1^  is  the  internal  energy  density  at  the  onset  of  melt 

IpDs  *s  the  «nerBy  required  to  complete  fusion  or  the  melting  process 
Iy  is  the  internal  energy  density  at  the  onset  of  vaporization 
IyAp  the  energy  required  to  complete  vaporization. 

As  the  diagram  indicates.  I  is  equal  to  1^  minus  internal  energy  lost  in  the 
fusion  and/or  vaporization  process.  Thus  pressure  due  solely  to  energy  addi¬ 
tion  is  not  allowed  to  increase  during  fusion  and  vaporization  until  these 
processes  are  complete.  Density  increases  will  still  provide  pressure  in¬ 
creases  through  the  first  term  in  the  equation. 

Metal  energy  at  the  onset  of  melt  increases  with  increasing  confine¬ 
ment.  Thus,  the  1^.  and  Iy  should  be  changed  as  pressure  increases.  HULL  uses 
the  GRAY  metal  (Reference  6)  fit  for  Ijj  as  a  function  of  volumetric  strain. 
The  function  is: 


Ijj  -  Iju  (1  +  aq  +  bqa) 


where : 
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a 


=  2  r0  -  2/3 
b  =  <r0  -  1/3)  (2  ro  +  1/3)  -  1 
n  =  i  -  pq/p 

Iw,,  =  internal  energy  density  at  onset  of  melt  at  ambient 
density 

We  increase  Iy  using  the  same  function,  i.e., 

I  =  — —  I 
V  *M0  V° 

Ipyg  and  Iy^p  ere  constant  regardless  of  the  zone's  compression. 

The  HULL  eqnati on-of-state  transitions  to  a  gamma-law  gas  at  internal 
energies  above  sublimation  energy  and  whenever  the  pressure  so  computed  is 
higher  than  the  Mie-Grunei sen  pressure.  The  pressure  comparison  ensures  a 
smooth  transition  between  equations  without  actually  determining  and  identi¬ 
fying  intersections  of  the  two  regions  in  P,  V,  I  space.  The  gamma-law 
equation  is: 


P  =  <Y  -  1)  P  (I  -  Is) 

where  Is  is  the  internal  energy  density  for  sublimation  and  y  is  Cp/Cy.  The  I 
in  this  equation  is  the  zone's  internal  energy — i.e.,  it  is  not  reduced  by 

IFUS  or  IVAP‘ 

There  are  a  number  of  constraints  applied  to  the  Mie-Grnneisen  pres¬ 
sure  to  more  physically  model  metal  behavior  ir  extreme  energy  and  density 
states.  The  constraints  are  more  easily  explained  by  expressing  the  Mie- 
Gruneisen  equation  as: 

P  =  Fy  (p)  +  F2  (p)  '  (I  -  I0>  +  P0 

When  material  in  a  zone  expands  under  tension,  it  should  not  be 
allowed  to  do  so  indefinitely.  Use  of  the  void  insertion  fracture  model  in 
HULL  can  prevent  extreme  expansion  states  from  occurring.  However,  this  model 
requires  dcta  which  is  often  unavailable.  In  these  cases  the  following 


72 


constraints  will  provide  some  physical  realism. 


The  product  ~  *o^  *s  l*m*te<*  if  H  is  less  than  -0.25.  The 

equation  employed  is: 

F2  (I  “  Io)  =  lF2  (I  “  Io)]  20  +  °-4>/3 

where  the  term  in  brackets  is  the  product  calculated  without  constraints.  The 
equation  provides  a  linear  descent  to  a  value  of  zero  for  the  product  at  ;  =  - 
0.4.  Because  of  an  IF  check,  the  equation  is  not  invoked  unless  ;  is  less 
than  -0.25.  The  point  of  controlling  this  product  is  simply  tbat  metals  at 
such  high  expansions  cannot  support  compression.  The  zone  contains  metal 
fragments,  and  such  fragments  cannot  support  compression  when  they  are  so 
highly  distended.  Without  this  control,  compressive  pressures  could  still 
exist  with  high  enough  I  values. 

The  Gruneisen  parameter  is  limited  so  as  not  1 3  exceed  2  r  .  In  highly 
expanded  states  the  Gruneisen  parameter  grows  without  limit  as  density  de¬ 
creases.  Since  the  Gruneisen  parameter  is  not  meant  to  be  applied  in  expanded 
zones,  it  is  reasonable  to  control  this  growth.  This  control  is  somewhat 
redundant  with  the  F^d— Iq)  control. 

If  the  calculated  Mie-Gruneisen  pressure  is  tensile,  there  are  several 
controls  which  are  applied. 

If  the  pressure  is  less  than  PHJN  (l-0.1p),  it  is  reset  to  PMIN  (1- 
O.lp).  Normally  P^j^  is  set  to  *  very  large  negative  number  (-1  x  10*°)  so 
that  this  control  is  not  invoked.  However,  the  same  equation  of  state  rou¬ 
tines  are  used  for  media  such  as  water  which  can  support  no  tension.  In  such 
cases,  this  control  is  required.  The  factor  (1-0, lp)  provides  a  variable 
pressure  as  p  changes.  The  product  Pjjjjq  (l~0.iy)  grows  slightly  more  negative 
as  p  decreases. 

If  the  Mie— Gruneisen  pressure  is  tensile  and  p  is  less  than  -0.2,  the 
pressure  is  reset  to: 


-20  (p  +  0.2)  (80  -  P)  +  P 


chec 


If  p  is  less  than  -0.25,  the 


where  P  is  the  pressnre  after  the  Pjjjn 
pressure  is  reset  to: 

100  (1  +  |J) 

These  two  controls  work  together  to  limit  total  pressnre  to  very  small  posi¬ 
tive  values  along  a  linear  ramp  from  P  =  0  at  p  =  -1  to  P  =  100  dynes/cm*  at  p 
=  0.  This  provides  essentially  zero  pressures  with  a  small  positive  slope  as 
density  increases. 

The  metal's  tensile  capability  is  limited  as  its  internal  energy 
density  approaches  melt  energy.  To  use  this  control,  the  quantity: 

IMT  =  +  °-5  ^US 

is  calculated.  This  is  also  the  quantity  used  to  limit  yield  strength  and 

will  fc:  discussed  later  in  this  section.  If  I  is  greater  than  0.6  Ijjy*  tlie 
calculated  tensile  Mie-Gruneisen  pressure  is  multiplied  by  the  factor: 


F  ~  <IjnrIo)(IinrI)  if  Ijfir  >  0,61my 
=0 

and  then  roset  to 

F  •  P  +  100  (1  +  JJ) 

This  control  limits  tensile  values  of  pressure  as  I  increases  to  I^y  amj 
phases  the  pressnre  into  the  distended  fragment-filled  zone  curve  previously 
discussed. 

Compressive  Mie-Gruneisen  pressures  are  limited  to  the  region  of  20  x 
10ia  dynes/cm*  (20  msgabars).  A  real  material  cannot  achieve  such  pressures 
without  also  having  internal  onergies  considerably  greater  than  I  jn  this 
case,  the  equation  of  state  will  have  switched  to  a  gamma-law  gas  which  has  no 
compressive  controls  on  it.  Therefore,  if  a  call  is  made  to  the  equation  of 
state  where  P  exceeds  20  megabars,  it  is  reset  to: 
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P  =  2  x  101*  +  2  X  10ai  (U-IK„)/(100  - 


where  pac  is  the  value  of  )J  when  P  is  2  x  101*.  This  provides  a  gentle  slope 
with  density  for  values  of  greater  than  )Ja0. 


Before  adjusting  yield  strength,  the  equation-of- state  routines 
provide  a  few  more  variables  required  in  other  sections  of  BULL. 


dP 

The  mixed-zone  equation-of- state  iterator  (Reference  5)  requires  (^ /  , 


,dP. 


where  x  is  the  specific  volume,  1/p.  HULL  calculates  (^)  and  multiplies  by 
-pa  to  obtain  this  quantity. 


The  adiabatic  bulk  sound  speed  is  required  for  several  purposes, 
including  time  step  control.  In  the  absence  of  heat  conduction,  the  adiabatic 
sound  speed  is  the  isentropic  sound  speed,  calculated  from  the  thermo¬ 
dynamic  identity  (Reference  1): 


_  ,ap.  _  ,dp.  p  „,ap, 
vs  op  S  op  I  p  91  p 


For  both  pressure  derivatives  HULL  takes  account  of  the  controls  invoked  and 
computes  derivatives  along  appropriate  paths. 

An  approximate  temperature  is  calculated,  for  edit  purposes  only,  from 
the  equatic  1: 


TCK)  -  I/C  +  T 
v  xo 

where  Cy  is  calculated  from: 

Cy  =  3R/AW 

R  =  gas  constant  =  8.134  x  107 
AV  =  the  metal's  atomic  weight 

and  where  To  i3  the  temperature  required  to  force  T  to  a  value  of  288°K  when  I 

ir.  I0* 

7. *5 
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B.  Yield  Streag tfe  Equtioai 

If  the  WORK  option  is  set  to  1  in  the  KEEL  (generator)  input  deck, 
yield  strength  is  varied  as  a  function  of  accumulated  plastic  strain  in  a  form 
as  seen  below. 


e 

P 

where  YLDST  is  the  initial  (zero  plastic  strain)  yield  strength  and  YLDMAX  is 
the  ultimate  yield  strength  which  occurs  when  ep  =  EPLAST.  If  WORK  =  0,  the 
ultimate  yield  strength  is  used  for  all  values  of  a 

The  calculation  of  ep,  incremental  strains  and  stresses,  and  resetting 
of  deviatoric  stresses  exceeding  the  yield  limit  are  carried  out  in  the  usual 
fashion  and  completely  described  in  the  previous  section.  Ihese  calculations 
will  not  be  repeated  here. 

After  an  ambient  energy  (possibly  work  hardened)  yield  strength  is 
calculated,  it  is  multiplied  by  a  thermal  softening  factor  which  varies  from  1 
at  ambient  energy  to  0  at  Ijjy.  Tjjy  is  the  melt  energy  at  which  all  yield 
strength  disappears.  HULL  sets  this  value  to: 

JMY  +  °*5IF0S 

where  1^  is  the  onset  of  aelt  energy  increased  as  a  result  of  increased 
density  and  Ipyg  the  energy  absorbed  during  the  melting  process.  This  is  a 
somewhat  arbitrary  assignment.  There  will  be  some  point  during  the  fusion 
process  when  enough  centers  of  melt  exist  to  drive  yield  strength  to  essen¬ 
tially  zero  values.  The  exact  position  of  this  point  is  unknown  and  will 
probably  never  be  known  since  observables  such  as  pressure  and  temperature  do 
not  vary.  It  does  not  seem  that  the  metal  would  have  to  be  entirely  melted  to 


essentially  zero  any  strength  characteristics.  On  the  other  hand,  when  a  few 
vary  small  centers  of  melt  are  being  created — i.e.,  at  the  onset  of  melt — it 
does  not  seem  that  the  strength  of  the  metal  would  be  appreciably  affected. 
To  avoid  both  extremes,  HULL  uses  the  midpoint  of  the  fusion  process. 


The  thermal  softening  factor  it  determined  from  a  trilinear  curve  as 


seen  below: 


THERMAL 

FACTOR 


where  the  pairs  (YFRAC1 , IFRAC1) ,  (YFRAC2 .IPRAC2 )  are  material  constants. 


HULL  employs  the  ratio  I/I^y  *nstead  of  temperature  to  calculate  a 
thermal  softening  factor  simply  because  tfco  relationship  between  temperature 
and  internal  energy  is  virtually  unknown  at  densities  other  than  ambient.  The 
HULL  technique,  in  effect,  implies  that  there  is  a  constant  relationship 
between  I/I^y  and  T/TM  at  anY  given  density  but  does  not  attempt  to  determine 


‘M 


the  value  of  the  constant. 


It  is  possible  that  a  highly  expanded  zone  could  have  a  significant 
value  of  yield  strength,  when,  in  reality,  the  rubble  in  the  zone  would  be  too 
distended  to  support  any  strength.  To  prevent  this,  yield  strength  is  multi¬ 
plied  by  another  factor  in  addition  to  work  hardening  and  thermal  softening. 
This  factor  is  1  for  U  values  greater  than  -0.2.  It  varies  linearly  from  1  at 
b  -  -0.2  to  0  at  b  =  -0.25  and  then  is  0  for  all  values  of  fJ  less  than  -0.25. 


HULL  does  not  calculate  an  explicit  multiplier  for  strain  rate  varia¬ 
tions  of  yield  strength.  Metals  demonstrate  the  property  of  having  an  almost 
constant  yield  strength  as  strain  rates  vary  from  10*  .ec”^  to  10s  sec  ^ . 
Since  conventioral  weapon  calculations  lie  almost  completely  within  this 
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ransc;  HDLL  assumes  that  the  strength  values  used  will  be  appropriate  within 
this  range. 

HULL  has  built-in  values  for  all  metal  material  properties.  Strength 
properties  can  be  changed  through  use  ot  the  MATPROP  option  in  the  HULL  input 
deck.  Other  property  changes  require  changing  DATA  statements  in  the  code. 
Tables  B-l  and  8-2  summarize  the  built-in  property  values  for  materials 
employing  the  Mie-<jruneisen  equation-of-state.  The  Hugoniot  data  was  obtained 
from  References  1,  10,  and  11  and  is  considered  very  accurate.  Melt  and 
vaporization  energy  values  were  obtained  from  many  sources  and  considered 
reasonably  accurate  unless  noted  otherwise  in  the  tables.  Yield  strength 
values  and  thermal  softening  data  were  also  obtained  from  many  sources 
including  material  property  tests  at  the  U.S.  Air  Force  Materials  Laboratory, 
the  fitting  of  cal cnl at ional  results  to  ambient  and  elevated  temperature 
cylinder  impact  tests,  and  the  fitting  of  calculatioDal  results  to  self 
forging  fragment  formation  experiments.  However,  users  should  consider  all  of 
the  strength  data  as  suspect  and  attempt  to  obtain  data  valid  for  the 
particular  metal  alloy  of  concern.  When  such  data  has  been  found  or 
generated,  the  user  can  override  built-in  strength  properties  by  entering  the 
following  cards  in  HULL  input:. 


MATFROP 


YLDST  =  Y1 
YFRAC1  =  YFj 
YFRA.C2  =  YF. 


ENDPROP 


YLDMAX  =  Y2 
IFRAC1  =  IF3 
IFRAC2  =  IF. 


EPLAST  =  E 


where  Y^,  Y^,  e,  YFj,  IF-,  YF2,  and  IFj  are  new  input  values,  and  N  is  the 
material  number  in  the  HULL  calculation.  The  example  illustrates  all  stress 
surface  values  being  changed.  Not  all  values  have  to  be  changed,  and  more 
than  one  material  can  be  altered. 


2.  Plain  Concrete  and  Sand  Equations  of  State 

HULL  contains  a  plain  concret’  equation-of-state  '.del  developed  to  pro¬ 
vide  reasonable  response  in  the  1-  to  300-kilobar  range.  It  is  based  on  a  very 
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TABLE  B-l.  HULL  HAIEUAL  FtOPKSTIES 


MATERIAL 

HULL 

NAME 

(gm?cc) 

c 

(10* cm/ sec) 

s 

r 

*  0 

PMIN 

(dynes/ cm*) 

Alum  in  tun 

AL 

2.71 

5.38 

1.337 

2.1 

-l.xlO20 

Ceramic 

(a12o3) 

CERAM 

3.9 

6.9 

1.45 

0.5 

-l.xlO1® 

Copper 

CU 

8.9 

3.958 

1.497 

2.0 

-l.xlO1® 

Iron 

FE 

7.86 

4.61 

1.73 

1.67 

-l.xlO*® 

Lead 

PB 

11.34 

2.092 

1  52 

2.0 

-l.xlO*® 

Nickel 

NI 

8.9 

4.8 

1.254 

1.83 

-1.X102® 

Nylon 

NYLON 

1.14 

2.29 

1.63 

0.87 

-10x10* 

Plastic 

(polyrubber) 

PLAST 

1.012 

0.854 

1.865 

0.9 

-10x10* 

RHA 

(thin  plate) 

RHA 

7.86 

4.61 

1.73 

1.67 

-l.xlO2® 

Stainless 

Steel 

SSTEEL 

SAME 

PROPERTIES 

AS  IRON 

Tantalum 

TA 

16.6 

3.423 

1.214 

1.69 

-l.xlO2® 

Teflon® 

TEFLON 

2.16 

1.34 

1.93 

0.9 

-10x10* 

Tungsten 

W 

18.1 

4.0 

1.268 

1.58 

-l.xlO20 

Uranium 

U 

19.05 

2.48 

1.53 

1.6 

-l.xlO20 

U-0.75K.TI 

UTI 

18.9 

2.48 

1.53 

1.6 

-l.xlO2® 

Water 

WATER 

1.0 

1.483 

1.75 

0.28 

0.0 

TABLE  B-l.  HULL  MATERIAL  PKOHSTIES  (CONCLUDED) 


JM0 

JFUS 

IV0 

IV0+IVAP 

y-1 (vapor 

MATERIAL 

(10*ergs/gm) 

(10*ergs/gm) 

(10*ergs/gm) 

(10*ergs/gm) 

(10»ergs/gm) 

equation) 

Allan  intun 

7.3 

3.96 

32.0 

100.0* 

119.2 

0.25 

Ceramic 

(A1203) 

9.55 

2.1 

50* 

100.0* 

100.0* 

0.25 

Copper 

4.78 

2.12 

16.3 

63.0 

53.1 

0.2  5 

Iron 

7.4 

2.74 

22.4 

86.8 

74.2 

0.25 

Lead 

0.67 

0.262 

2.96 

11.2 

9.4 

0.25 

Nickel 

6.6 

3.1 

20.1 

84.7 

72.9 

0.25 

Nylon 

1.9 

2.1 

50.0* 

100.0* 

100.0* 

0.25 

Plastic  20.0* 

(polyrubber) 

2.0* 

50.0* 

100.0* 

100.0* 

0.25 

RHA  7.4 

(thin  plate) 

2.74 

22.4 

86.8 

74.2 

0.25 

Stainless 

Steel 

SAME 

PROPERTIES  AS 

IRON 

Tantalum 

4.32 

1.59 

12.5 

48  3 

43.2 

0.25 

Teflon® 

20.0* 

2.1* 

50.0* 

100.0* 

100.0* 

0.25 

Tungsten 

4.77 

1.84 

13.6 

58.4 

46.0 

0.2 

Uranium 

1.4 

0.49 

5.5 

28.2 

22.7 

0.25 

U-0.75%TI 

1.4 

0.49 

5.5 

28.2 

22.7 

0.25 

Water 

7.1 

3.35 

23.8 

44.7 

44.7 

0.25 

•NO  DATA, 

ARBITRARY  VALUES  ASSIGNED 

TABLE  B-2.  HULL  STSB1GTB  PROPERTIES 


MATERIAL 

YLDST 

(10* 

DYNES/ CM*) 

WYLDST 

(10* 

DYNES/ CM*) 

EPLAST 

YF1 

IF1 

IF2 

EF2 

SHEAR 

MODULUS 

(10** 

DYNES/CM*) 

Aluminum 

2.9 

2.9 

0.0 

0.9 

0.5 

0.9 

0.5 

2.69 

Ceramic 

(AlaO,) 

80.0 

80.0 

0.0 

0.21 

0.45 

0.05 

0.75 

10.0 

Copper 

1.2 

4.5 

0.3 

0.78 

0.41 

0.78 

0.41 

4.64 

Iron 

4.69 

5.5 

0.3 

0.9 

0.5 

0.9 

0.5 

6.41 

Lead 

0.3 

0.3 

0.0 

0.9 

0.8 

0.9 

0.8 

1.113 

Nickel 

6.0 

6.0 

0.0 

0.8 

0.5 

0.8 

0.5 

8.9 

Nylon 

0.5 

0.5 

0.0 

0.9 

0.5 

0.9 

0.5 

0.368 

Plastic 

(Polyrubber) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

RHA 

15.0 

15.0 

0.0 

0.9 

0.5 

0.9 

0.5 

6.41 

Stainless 

Steel 

6.89 

10.0 

0.3 

0.9 

0.5 

0.9 

0.5 

7.3 

Tantalum 

5 

5 

0.0 

0.5 

0.4 

0.5 

0.4 

6.56 

Teflon® 

0.5 

0.5 

0.0 

0.9 

0.5 

0.9 

0.5 

0.233 

Tungsten 

20.0 

20.0 

0.0 

0.9 

0.5 

0.9 

0.5 

16.0 

Uranium 

8 

15 

0.5 

0.9 

0.5 

0.9 

0.5 

5.1 

U-0.75%TI 

8 

15 

0.5 

0.9 

0.5 

3 

0.5 

5.1 

Water 

0.0 

0.0 

0.0 

0.0 

O.t 

0.0 

0.0 

small  amount  of  data,  since  most  available  data  is  at  stress  levels  well  below 
1  kilobar.  However,  the  equation-of- state  has  been  used  successfully  in 
calculations  of  penetrations  into  concrete  at  velocities  from  a  few  hundred  to 
over  1,000  feet/second  (Reference  12). 

The  model  was  constructed  primarily  from  Hugoniot  data  on  one  specific 
concrete  (Reference  13)  and  static  yield  strength  data  on  several  concretes 
(Reference  14).  Assumptions  were  made  for  a  value  for  Poisson's  ratio,  a 
shear  modulus  formula,  a  plastic  flow  rule,  and  other  parameters. 

Gregson  (Reference  13)  generated  Hugoniot  data  to  stress  levels  over  500 
kilobars  for  a  granite  aggregate  concrete.  The  aggregate  had  an  average 
diameter  of  1/8  inch.  The  initial  density  varied  from  mix  to  mix  but  averaged 
2.185  gm/cc.  The  solid  density  had  an  average  of  2.673  gm/cc.  The  concrete 
consisted  of  18  percent  voids,  25  percent  granite  aggregate,  20  to  25  percent 
qnartz  grains,  and  25  to  30  percent  cement  paste  by  volume.  The  small  aggre¬ 
gate  size  was  dictated  by  his  shock-generating  method,  which  consisted  of  gas- 
gun  launched  flyer  plates.  He  employed  several  measurement  techniques  to 
establish  stress  and  velocity.  The  final  data  was  converted  to  stress  versus 
excess  compression,  p. 

The  lowest  point  in  °reg«on's  data  is  at  0.75  kilobar.  This  is  a  fairly 
consistent  precursor  value  seen  in  his  stress  versus  time  data.  BULL  uses 
that  point  to  define  the  density  ratio  at  which  significant  cracking  and  void 
filling  occurs. 

The  estimated  unloading  path  is  simply  one  which  allows  the  concrete  to 
return  to  its  solid  density  when  loaded  past  the  point  where  all  voids  are 
filled.  It  is  a  straight  line  which  connects  with  the  loading  curve  at  the 
total  void  filled  point  with  a  slope  equal  to  the  loading  curve  at  that  point. 
There  is  no  data  for  this  line  beyond  the  solid  density  point  at  zero  stress. 

The  yield  strength  of  concrete  increases  as  confining  pressure  is  in¬ 
creased,  The  only  data  available  at  high  confining  pressures  is  static  data 
generated  by  Chinn  (Reference  14).  Chinn's  results  were  used  along  with  some 


typical  low  stress  level  results  from  McHenry  (Reference  15).  The  yield 

strength  and  pressure  are  normalized  by  f',  the  concrete's  unconfined 

compressive  strength.  Two  theoretical  points  were  chosen.  One  defines 

failure  (flow)  at  the  unconfined  compressive  strength  point,  i.e.,  Y  =  at  P 

-  3./3  f*  The  other  is  included  because  a  minimum  pressure  criterion  is  used 

to  define  failure.  This  minimum  pressure  was  chosen  as  that  occurring  in  an 

unconfined  tensile  test  assuming  failure  when  stress  is  0.1  f*  (Reference  16). 

c 

It  is  numerically  pleasing  to  set  Y  =  0  at  this  failure  point. 

The  yield  strength  curves  are  approximated  in  HULL  with  three  straight 
line  segments: 

Y  =  0  if  P  <  -O.lf* 

—  c 

Y  =  3  (P  +  O.lf * /3)/l .1  if  -O.lf * /3  <Pi  f'/3 

c  c  —  c 

Y  =  P  +  2/3  f*  if  f'/3  <  P  <  30f * 

c  c  c 

Y  =  30.67  f*  if  P  >  30  f' 

c  c 

the  limit  at  P  =  30T  is  not  indicated  in  Chinn's  data.  It  is  included 
because  it  is  felt  there  should  be  a  limiting  value.  It  is  included  at  the 
level  shown  simply  because  it  is  not  known  to  exist  at  a  lower  level. 

It  should  be  emphasized  that  the  yield  strength  data  was  developed  in 
static  tests.  Higher  yield  strengths  may  exist  at  strain  rates  of  interest  in 
conventional  munitions  work,  but  data  simply  does  not  exist. 

Using  the  previous  equations  for  Y,  Gregson’s  Hugoniot  stress  data  was 

reduced  to  provide  an  Hugoniot  pressure  versus  excess  compression  curve  by 

simply  subtracting  2/3  Y  from  the  stress  data  assuming  an  f'  value  of  5,000 

c 

PSI. 

Below  the  precursor  pressure  level  of  0.358  kilobar  (1/3  x  0.75  kilobar), 
Hugoniot  pressure  is  fit  with  a  straight  line: 

PH  =  V 

where  Kq  is  determined  from  the  empirical  Young's  modulus  relationship 
(Reference  10) 


83 


i 


E  =  33I1*5  (fJ.)°‘5PSI  (W  in  lb/ft*  =  density) 
and  the  relationship 

K  =  E/3(l-2or) 

with  a  value  of  Poisson's  ratio,  a,  chosen  as  0.2. 


Designating  the  value  of  p  at  0.3S8  kilobar  as  Pg,  the  loading  Hugoniot 
pressure  was  fit  with  the  following  equations: 


PH(Kb) 


RANGE  OF  p 


V 

0.358  +  78.62  (p-p£) 

8  +  130  (p  -  0.1) 

21  +  420  (p  -  0.2) 

63  +  650  (p  -  0.3) 

76  +  625.5  (p  -  0.32) 

+  1720.19  (p  -  0.32) 

150  +  39.68  (p  -  0.414) 

166  +  360.5  (p  -  0.572) 
+  864  (p  -  0.572) 


-0  to  p£ 

PE  t0  o*1 
0.1  to  0.2 
0.2  to  0.3 

0.3  to  0.366  (lockup  point) 


0.366  to  0.414 

0.414  to  0.54  (solid 

-  solid  phase  transition) 

0.54  to  ® 


The  maximum  density  unloading  curve  was  fit  by  the  equation: 

PH  =  784  (p  -  0.223) 

HULL  maintains  a  zone  variable,  Pg,  for  concrete  which  defines  the  maximum 
value  of  p  previously  seen  by  the  zone.  For  values  of  Pg  less  than  P£* 
loading,  unloading  and  reloading  ell  occur  along  the  same  elastic  path.  For 
values  of  Pg  greater  than  p  at  lockup  (0.366)  loading,  unloading,  and 
reloading  take  place  along  the  same  path.  For  Pg  values  between  these  ex¬ 
tremes,  unloading  and  reloading  up  to  the  virgin  loading  curve  take  place 
along  lines  where  K  is  chosen  to  vary  linearly  between  KQ  and  784  kilobars. 
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The  shear  modulus  is  calculated  assuming  that  Hooke's  lav  holds.  That  is: 


G  = 


3(l-2o) 

2(l+o) 


* 


where  is  the  bulk  modulus,  p  (^) ,  and  o  is  Poisson's  ratio. 

Deviators  are  reset  as  though  the  yield  surface  were  a  von  Mises  surface — 
i.e.,  not  a  function  of  pressure.  The  resulting  plastic  flow  rule  is  incon¬ 
sistent  with  the  Mohr-Coulomb  yield  surface  proposed  for  the  concrete.  How¬ 
ever,  the  alternative,  such  as  a  moving  CAP  model  (Reference  2),  is  computa¬ 
tionally  complex  and  cannot  be  justified  considering  the  source  of  the  yield 
strength  model.  Strain  rate  effects  are  probably  far  more  important  than  a 
consistent  flow  rule. 

The  Hugoniot  pressure  previously  discussed  could  be  combined  with  a 
Gruneisen  parameter  into  a  complete  equation  of  state.  However,  Gruneisen 
parameter  measurements  for  concrete  are  typically  around  0.1.  For  this  rea¬ 
son,  the  equation-of-state  ignores  the  Gruneisen  parameter  and  Hugoniot  pres¬ 
sure  is  used  as  total  pressure.  This,  of  course,  means  that  internal  energy 
changes  have  no  effect  on  pressure,  and  the  equation-of-state  could  not  be 
used  for  energy  deposition  problems. 

The  variable,  f£,  is  set  to  the  CGS  equivalent  of  5,000  PSI.  If  other 
strengths  are  desired,  it  can  be  changed  in  the  data  statement  that  now  sets 
it.  Changing  f£  will  vary  the  fracture  point  and  the  yield  surface.  It  will 
have  no  effect  on  the  pressure. 

The  sand  equation-of-state  was  developed  to  simulate  a  dry  40  percent 
porosity  quartz  sand.  It  combines  low  pressure  hydrostatic  data  from  Refer¬ 
ence  14  with  pure  quartz  Hugoniot  data  from  Van  Thiel  (Reference  10).  The 
sand  model  assumes  that  strength  is  negligible  at  pressure  levels  of  interest. 
This  is  a  valid  assumption  as  long  as  strength  levels  of  a  few  hundred  PSI  are 


unimportant  in  a  calculation.  There  is  evidence  that  sand  can  achieve  such 
small  levels  o£  strength  when  highly  crushed. 

Sand  has  a  very  small  Gruneisen  parameter  which  can  be  ignored.  The 
equation-of-state,  then,  is  simply  a  pressure  versus  excess  compression  fit. 
The  virgin  loading  curve  is  described  by  a  very  small  elastic  segment: 

P  =  K0H 

from  (i  =  0  to  p  =  0.0084.  KQ  is  computed  from  the  initial  density  of  1.6 
gm/cc  and  the  initial  sound  speed  of  0.61  x  10*cm/sec.  The  pressure  level  at 
the  end  of  the  elastic  point  is  then  0.05  kilobar  (approximately  50  atmos¬ 
pheres)  the  elastic  segment  is  followed  by  two  linear  crushing  regions: 

P  =  0.05  +  4 . 9  6  JJ  (Kilobars) 

for  pi  from  0.0084  to  0.2 
and  P  =  1.0  +  73.3pJ  (Kilobars) 

for  pi  from  0.2  to  0.425 

At  a  pj  of  0.425  and  a  pressure  of  17.5  kilobars,  the  sand  is  assumed  to  have 
crushed  out  all  voids  and  the  loading  curve  follows  the  quartz  Hugoniot:. 

P  =  2.98q/ (l-2.33q) *  (Kilobars) 

where : 

V  =  1  “  Pc/P 

To  avoid  the  pole  in  this  Hugoniot,  the  Hugoniot  continues  as  a  constant 
modulus  curve  at  q  =  0.98/s. 

Unloading  and  reloading  below  the  lockup  point  of  17.5  kilobars  and  above 
the  elastic  point  of  0.05  kilobar  follows  a  straight  line  curve  which  has  a 
bulk  modulus  continuously  varying  from  that  at  the  lockup  point  to  that  along 


The  sand  is  not  allowed  to  sustain  any  tension. 

Both  concrete  and  sand  require  an  extra  zone  variable,  the  maximum  value 
of  p  ever  seen  by  the  material.  To  successfully  use  these  models,  the  user 
must  set  NHIST  =  1  in  the  KEEL  input  deck.  This  will  cause  PLANK  to  increment 
NH1ST  by  one  which  will  allot  one  extra  variable  per  zone  after  the  variables 
required  for  carrying  stresses.  This  will  hold  *or  concrete  and/or  sand. 

3,  Miscellaneous  Solid/Liquid  Equationr-of— State 

There  are  a  few  equations-of-state  in  HULL  other  than  those  following  the 
Nie-Gruneisen  formulation  (the  materials  in  Tables  B-l  and  B-2)  and  the  con¬ 
crete  and  sand  models.  These  will  now  be  briefly  discussed. 

Undetonated  AFX-108  is  very  simply  modelled  for  use  in  low  velocity  explo¬ 
sive-filled  bomb  and  warhead  impacts.  The  material  name  in  HULL  is  AFX1.  The 
Hugoniot  pressure  is  calculated  from: 

PH  "  »1  *  *« 

where  js  set  to  10.3  x  10*  dynes/cma.  The  Gruneisen  parameter  is  set  to  a 
constant  value  of  1J.  and  pressure  is  computed  from  the  Mie-Gruneisen  equation 
with  no  controls  or  tensile  values.  The  only  control  is  a  compressive  one.  )J 
is  not  allowed  to  exceed  0.9.  This  strange  equation^of-state  is  the  result  of 
having  only  a  single  loading  modulus  for  the  material. 

Tillotson  equations-of-state  exist  in  HULL  for  granite  and  tuff.  They 
were  developed  for  nuclear  cratering  calculations  and  should  not  be  considered 
valid  at  low  conventional  weapons  research  pressure  levels.  They  have  not 
been  exercised  with  FLUXES  =  3,  and  the  user  is  advised  to  avoid  them  if 
possible. 

4.  Air  Equations  of  State 

Two  air  equations-of-state  are  implemented  in  EOS  6.  If  t>e  user  simply 
asks  for  AIR,  the  equation-of- state  uses  the  strange  form: 

P  -  (Y-l>  >IP0I0 
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where  y  -  1.4  and  Io  is  set  to  return  one  atmosphere  pressure  at  p  =  pQ.  The 
purpose  of  this  form  is  to  j  *ovide  air  at  approximately  the  correct  density 
but  with  a  very  low  sound  speed.  Small  amounts  of  air  are  occasionally 
trapped  in  zones  as  the  result  of  motion  of  more  substantial  media.  High 
pressures  in  this  trapped  air  can  drive  the  time  step  to  very  small  values  and 
dominate  the  problem,  when  the  pressure  in  the  air,  and  perhaps  even  the 
presence  of  the  air  itself,  is  unimportant. 

If  correct  air  pressures  are  important,  then  the  user  must  set  AIREOS^T  as 
a  SAIL  option.  AIRE0S=1  insures  that  the  DOAN-NICKEL  variable  gamma  air  model 
will  be  used.  This  is  a  highly  accurate  model  over  broad  ranges  of  density 
and  energy  and  is  described  in  detail  in  Reference  26. 

Because  the  code  was  used  to  simulate  nuclear  fireball  calculations  at 
various  heights  above  sea  level,  HULL  is  capable  of  simulating  several  atmos¬ 
pheres  (air  initial  conditions).  These  are  described  in  Reference  27.  If 
ATMOS  is  not  specified,  a  constant,  sea  level  atmosphere  will  be  assumed. 

5.  High  Explosive  Equation  of  State 

There  are  several  high  explosive  equations  of  state  used  in  HULL.  The 
detonation  products  of  standard  military  explosives  are,  for  the  most  part, 
represented  by  the  JWL  formulation  (Reference  18).  In  this  formulation, 
pressure  is  computed  from: 

P  «  A  (1  -  £3)  e  -(Rl/q)  +  B  (1  -  £3)  e  -<R2'q>  +  Upl 

K1  R2 

where  q  =  p/pQ,  I  is  internal  energy  per  unit  mass  and  A,  B,  Rt,  R*  and  u>  are 
material  constants. 

This  equation-of-state  is  implemented  for  Composition  B,  Octol  (Reference 
22)  and  C4  (Reference  23).  JWL  data  exists  for  implementation  for  several 
more  explosives  (References  18-23), 

The  JWL  equation  has  proven  especially  useful  in  accurately  predicting 
metal  acceleration.  Its  utility  largely  derives  from  the  manner  in  which  the 
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material  constants  are  determined.  They  are  forced  to  fit  several  consistency 
criteria  snch  as  total  deliverable  energy  equal  to  that  measured  calorimetri- 
caliy,  pressure  at  Chapman-Jouget  (CJ)  density  equal  to  measured  CJ  pressures, 
etc.  The  remaining  parameters  (Ra,  Ra>  and  «)  are  fit  through  an  iterative 
experimental/calculational  procedure  involving  measured  velocities  of  cylin¬ 
ders  and  spheres. 

At  low  densities,  the  last  term  in  the  JWL  equation  dominates.  HULL  takes 
advantage  of  the  calculational  simplification  offered  by  this  singlo  term  at 
high  expansions.  Whenever  density  of  the  detonation  products  decreases  to 
one-tenth  initial  density,  pressure  is  computed  from  the  single  pul  term. 

For  densities  greater  than  one- tenth  ambient,  the  JWL  pressure  does  not 
proceed  to  zero  as  energy  proceeds  to  zero.  This  situation  is  avoided  in  HULL 
by  decreasing  pressure  whenever  energy  is  less  than  1  x  10*  ergs/gm.  The 
calculated  pressure  is  multiplied  by  1/1  x  10*.  This  arbitrary  factor  allows 
low  energy,  relatively  high  density  gases  to  react  more  realistically.  By  the 
time  energy  has  dropped  from  initial  4  to  5  x  10x0  erg/gm  levels  to  1  x  10», 
the  explosive  has  completely  delivered  its  energy  to  snrrounding  media  and  is 
merely  expanding  to  fill  in  the  available  grid  regions. 

The  JWL  equation-of-state  calculates  a  temperature  for  editing  purposes 
only  from  the  equation: 

T(°K)  =^1+^ 

where  i8  1.65  x  10  ®  and  is  1375.  These  constants  were  found  to  fit 
available  data  for  Composition  B  and  are  used  for  all  JWL  explosives. 

The  JWL  constants  A,  B,  R^,  and  o>  are  designated  as  Bj,  B2,  B3,  B4,  and 
Bj,  respectively,  in  HULL.  Values  of  these  constants  are  listed  in  Table  B-3. 

HULL  can  provide  burn  progression  (in  2D  only)  based  on  the  CJ  condition 
(D=C+U)  and  the  direction  of  pressure  gradients  in  the  unburned  explosive 
(BURN=1)  or  a  programmed  burn  can  be  used  (BURN=2).  A  programmed  burn  is 
suggested  when  an  accurate  detonation  velocity  is  desired  and  when  zone  sizes 
are  such  that  pressure  at  the  burn  front  is  significantly  less  than  Cj  pres¬ 
sure.  Values  of  detonation  velocity,  D,  are  included  in  Table  B-3  for  use  in 
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BURN-2.  Pressure  and  y^j  =  Cp/Cy  values  at  the  CJ  front  are  also  included  in 
the  table  for  calculational  comparison  purposes. 

If  BURN=2  is  desired,  the  user  inserts  values  for  TDET  (time  that  the 
explosive  region  begins  burning),  VDET  (detonation  front  velocity),  and  XDET, 
YDET,  and  ZDET  (ignition  point).  These  values  are  inserted  as  the  last  ele¬ 
ment  in  the  PACKAGE  setting  up  the  explosive  region. 


TABLE  9-3 .  JWL  EXPLOSIVE  CONSTANTS  (OGS  UNITS) 


BULL 


EXPLOSIVE 

NAME 

p0 

J0 

A 

B 

«1 

%  “ 

COMP  B 

(64  RDX,  36  TNT) 

CBBRN 

1.72 

4.92  x 
10“ 

5.242  x 
10“ 

0.07678  x 
10“ 

4.2 

1.1  0.34 

OCTOL 

OCTBRN 

1.82 

5.27  x 
10“ 

7.486  x 
10“ 

0.1338  x 
10“ 

4.5 

1.2  0.375 

C4  C4BRN 

(91  RDX,  5.3 
DI(Z-ETHYLHEXYL) , 

2.1  POLYISOBUTYLENE, 

1 .6  MOTOR  OIL) . 

1.601 

5.62  x 
10“ 

6.0977  x 
10“ 

0.1295  x 
10“ 

4.5 

1.4  0.25 

EXPLOSIVE 

NAME 

D 

C 

ycj 

P  . 

CJ 

COMP  B 

CBBRN 

7.98 

x  105 

0.01082  x 

10“  2.706 

295 

x  10» 

OCTOL 

OCTBRN 

8.48 

X  10* 

0.01167  X 

10“  2.83 

342 

X  10» 

C4 

C4  BRN 

8.193 

x  10* 

0.01043  x 

10“  2.838 

280 

x  10» 

It  is  sometimes  convenient  to  calculate  the  work  performed  by  an  explosive 
in  expanding  from  the  CJ  point  to  some  density  p. 

The  JWL  adiabat  (isentrope)  is  given  by: 
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P  =  A  c  +  B  e  +  C(n)“+1 


the  values  of  C  for  this  adiabat  are  alsc  given  in  the  table.  Integration  of 
-P  dp/p1  from  p£j  to  the  density  of  interest  will  provide  deliverable  energy. 

Other  detonation  product  equations-of-sta te  exist  in  HULL.  They  can  be 
nsed  (cautiously)  for  ANFO,  BARATOL,  TNT,  PBX9404,  and  Pentolite.  They  use 
various  forms  for  the  equation-cf-state  and  were  developed  for  nuclear  explo¬ 
sive  simulation  purposes. 

ANFO,  TNT,  and  Pentolite  detonation  products  pressures  were  fit  to  the 
LSZK  equation  (References  24  and  25) 

P  =  BlPi  +  B2pB3 

with  the  parameter  values  as  seen  in  Table  B-3.  The  TNT  and  Pentolite  data 
are  considered  to  be  very  accurate  at  large  expansions  and  were  extensively 
tested.  The  equations  have  not  been  checked  for  near  Cj  behavior  accuracy 
which  would  be  dominant  in  metal  acceleration  problems.  We  suggest  applica¬ 
tion  of  the  JWL  equation  for  this  region.  The  ANFO  data  used  in  HULL  was  from 
a  preliminary  ANFO  model  developed  by  Lawrence  Livermore  Laboratory.  Better 
ANFO  data  now  exists  and  is  contained  in  Reference  20.  HULL  has  not  been 
updated  because  of  the  lack  of  interest  in  ANFO  in  conventional  munitions 
work. 


The  Baratol  equation-of-state  uses  a  gamma  law  formulation  simply  because 
better  data  does  not  seem  to  be  available.  The  material  parameters  for  this 
explosive  are  also  given  in  Table  B-3.  Based  on  comparisons  with  other  explo¬ 
sives,  we  can  state  that  the  constant  gamma  assumption  will  provide  metal 
velocities  close  to  20  percent  too  high. 

A  PBX9404  fit  to  the  JWL  equation-of-state  is  included  in  HULL.  The  fit 
was  further  modified  to  reproduce  airblast  data  in  the  range  of  15  kilobars 
and  below.  It  should  not  be  used  for  metal  acceleration. 

The  predetonation  equatioas-of-state  in  HULL  mostly  use  the  form: 
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P  =  Bjji  +  BjJ*  +  B3;»  +  apl 

Tbit  provides  a  simple  approx last ion  to  pressure  for  solids  expected  to  load 
along  a  Hugoniot  and  for  which  subsequent  unloading  and  reloading  are  on  im¬ 
portant.  The  pressure  is  not  allowed  to  become  negative  in  this  equation-of- 
state. 

Temperature  is  calculated  fox  edit  purposes  only  and  is  fit  using  the 
equation: 

T  =  ClI  +  if  T>288#K 
T  =  C3I  if  T<288®£ 

Values  of  the  coefficients  used  in  the  pressure  and  temperature  calcula¬ 
tions  are  included  in  Table  B-4. 

TABLE  B-4.  PIEQBTORAXHM  EXPLOSIVE  PA1AMBTEKS  (CBS  UNITS) 

HULL 


EXPLOSIVE 

NAME 

B1 

®2 

®3 

Cl 

C2 

C3 

COMP  B 

COMP  B 

1.282 

X 

0 

1.1932  x 

3.36  x 

200 

1.101  x 

10** 

10** 

10~8 

10 

OCTOL 

OCTOL 

1.282 

X 

0 

1.1932  x 

3.36  x 

200 

1.101  X 

10** 

10** 

10  8 

10  7 

TNT 

TNT 

1.282 

X 

0 

1.1932  x 

3.312  x 

200 

1.085  x 

10** 

10'* 

10~8 

10_/ 

PENTOLITE 

PEN 

1.319 

X 

3.614  x 

2.295  x 

3.312  x 

200 

1.085  x 

10** 

10** 

10** 

10-8 

10  7 

ANFO 

ANFO 

2.4 

X 

1.08  x 

7.2  x 

3.3  x 

200 

1.081  x 

io*« 

10** 

10*® 

10"8 

10  7 

The  Baratol  (HULL  nave  BARTOL)  predetonation  equation  is  a  Mie-Gruneisen 
form.  The  coefficients  are: 

CQ  =  2.4  x  10* 
s  =  1.66 


o  =  °-2 

Temperature  fits  the  fora: 

T  -  Cjl  +  Cj 

where  -  3.36  x  10  ®  and  C2  =  200. 

6.  High  Pressure  Geologic  Equatioir-of-State 

A  high  pressure  equation-of-state  for  geologic  materials  was  added  to  BULL 
to  provide  the  necessary  material  properties  for  the  calculation  of  meteor 
craters.  It  is  designed  primarily  to  provide  a  good  fit  to  Hugoniot  data 
between  50  kilobars  and  1  megabar.  Many  details  of  the  Hugoniot  are  only 
grossly  approximated  below  50  kb  and  so  the  equation-of-state  should  be  used 
with  caution  in  this  area.  It  is  based  on  the  information  in  References  28 
and  29  modified  to  aliow  a  fit  to  a  Hie— Gruneisen  equation-of-state. 

The  Hugoniot  is  modelled  in  several  connected  sections.  The  Hugoniot 
lock-up  curve  is  the  curve  which  describes  material  behavior  for  a  geologic 
material  with  no  voids.  The  crush-up  curve  describes  how  a  material  with  a 
non-zero  void  fraction  behaves  prior  to  reaching  a  zero-void  condition.  In 
this  formation.  is  the  density  at  ambient  conditions  and  pr  is  the  density 
at  ambient  pressure  for  the  zero-void  material.  The  crush-up  loading  curve  is 
described  by  a  single  bulk  modulus.  K^.  The  lock-up  loading  curve  is 
described  in  several  segments.  A  constant  bulk  modulus  segment  connects  the 
1-atmosphere  point  with  the  point  at  which  the  crush-up  curve  meets  the  lock¬ 
up  curve  (jj^,  P^).  This  constant  bulk  modulus  curve  then  turns  into  a 
variable  modulus  curve  between  and  Pgj.  The  curve  is  described  by  the 
Schuster-Isenberg  form 

PH  =  Vr  -  ( VKo>  (l-e-^*) 


93 


[  where  Pr  =  P/Pj.  -  i  and  KQ  end  are  initial  and  maximum  bulk  moduli  and 

|  is  a  parameter  which  determines  the  density  range  between  the  points  at  which 

|  and  Km  dominate.  At  the  point  (Pgp  P|jt)  tlie  Hugonlot  shifts  again  to  a 

i  constant  modulus  form 

t 

PH  =  *HI  (Mr  "  >JHI)  +  PHI 

the  unloading  path  taken  depends  upon  u  which  is  the  maximum  value  of  p 

UftZ 

seen  by  the  material.  If  Pfflax  is  greater  than  then  all  unloading  and 

reloading  occurs  along  the  single  lock-up  curve.  If  p  is  less  than  p, , 

max  l 

then  unloading  and  reloading  occur  along  a  path  parallel  to  the  lock-up  curve. 
The  Hugoniot  pressure  is  unlimited  in  tensile  extent,  but  the  total  pressure 
is  limited  to  values  greater  than  Pajfl»  where  Pfflin  *s  a  material  parameter. 
In  addition,  total  tensile  pressure  is  limited  as  specific  internal  energy,  I, 
approaches  I  »  the  value  at  which  strength  effects  are  no  longer  important. 
In  this  case  the  total  tensile  pressure  is  multiplied  by 

2.5  Max  (I  -  Max  (I,  0.6T  ),0)/I„ 

ffi  m  ID 

This  formulation  insures  that  tensile  pressure  values  which  have  passed  the 
P_,_  test  are  not  modified  unless  I  exceeds  0.6  I_. 

OlO  m 

The  total  pressure  is  calculated  using  the  Mie-Gruneisen  form 

P  =  PH  (i  -  Fp/2)  +  Tpl  +  p0 
where  F,  the  Gruneisen  parameter,  is  given  by 

r  =  WP 


subject  to  the  restriction  that  F  cannot  exceed  2  T0.  PQ  i*  given  by 

„  1.013X10*  ■  IT  T  \  A  /I 

p  =  - - -  mm  (i,  io)  dynes/cm* 

o  lo 

where  Io  is  ambient  (1  atmosphere)  specific  internal  energy.  PQ  then  insures 
that  material  pressure  will  initially  be  1  atmostphere.  The  total  solid  or 
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liquid  pressure  P,  is  limited  by  P«IN  and  an  Ia  factor  as  discussed  previous¬ 
ly.  In  addition,  it  is  limited  to  very  small  positive  values  if  JJ  is  less 
than  -0.25  (the  theory  being  that  it  is  essentially  rubble  unable  to  support 
tension  at  these  low  densities). 

In  addition  to  P,  a  gas  pressure,  Pg,  is  calculated  if  I  exceeds  lm.  This 
pressure  is  calculated  using  a  gamma  lav  formulation 

p6  =  (Y-i)  P  d-y 

PG  is  compared  to  the  solid/liquid  pressure,  P,  and  P  is  reset  to  the  maximum 
of  the  two  values. 

The  shear  modulus,  G,  is  calculated  from 

G  =  (G  u  -  (G  -G  )  (1-e^G) )  (l-I/I  ) 

at  no  in 

if  1  <  I 

m 

and 

G  =  0  if  I  2  I 

m 

where  Gc,  Gm,  and  jiG  are  material  parameters. 

The  yield  strength,  Y,  is  calculated  from 

Y  =  Yr* Yj'Min  (A  +  BP,  Yasx> 
where  Yf  is  a  rubble  factor  set  as  follows 

Y  =  1  if  JJ  >  -0.2 

r 

Tr  =  0  if  JJ  <  -0.25 

Yr  =  20  (»J  +0.25)  if  -0.25  <  JJ  <  -0.2 

The  factor  insures  that  rubble  cannot  sustain  any  yield  strength.  Yj.  is  an 
internal  energy  factor  and  is  given  by 

YI  -  0  if  1  2  I„ 

Yj  =  1  -  I/IB  if  I  < 
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This  is  a  thermal  softening  factor.  The  term  A  +  BP  provides  a  Mohr-Coulomb 


yield  strength  which  reaches  saturation  at  Y  =  Yaax* 

Temperature  is  calculated  for  edit  purposes  by  arbitrarily  setting  ambient 
internal  energy  to  1  x  10*  ergs/gm  and  computing  a  Cy  value  to  return  288  K  at 
ambient  energy. 

The  attached  Table  B-5  lists  materials  which  utilize  this  equation-of- 
state,  the  HULL  names  for  the  materials,  and  the  material  quantities  built 
into  HULL  for  the  materials. 


TAM*  B-5.  HIGH  PHES STOLE  GEOLOGIC  PtOPKHTIES 


Dry 

Dry 

Dry 

Saturated 

Moenkopi 

Kaibab 

Coconino 

Coconino 

MATERIAL 

Sandstone 

Limestone 

Sandstone 

Sandstone 

Limestone 

HULL  NAME 

MOENK 

KAIBAB 

DRYCOC 

WETCOC 

LIMEST 

po(gm/cc) 

2.25 

2.30 

2,03 

2.35 

2.58 

Pr(gm/cc) 

2.68 

2.70 

2.68 

2.35 

2.70 

Kj( dynes /cm*) 

60.  E9 

100. E9 

60.  E9 

225. E9 

500. E9 

K0(dynes/cma) 

450. E9 

600. E9 

450. E9 

225. E9 

600. E9 

dynes /cm*) 

700. E9 

800. E9 

700. E9 

580. E9 

800. E9 

KpT(dynes/cm*) 

4.656E12 

2.942E12 

4.656E12 

2.22E12 

2.25E12 

ui 

^Lock 

> 

0.029208 

0.03541 

0.0488 

0 

0.196 

0.65 

0.2963 

0.7:>373 

0.74468 

0.33333 

0.25 

0.4 

0.25 

0.3 

0.4 

PG 

0.0015 

1.0 

0.0015 

0.0015 

1.0 

p  ,  (dynes/cm*) 

j/*tn 

(Synes/cm*) 

-0.667E8 

-0.15E9 

-0.667E8 

-0.667E8 

-0.15E9 

13.55E9 

22.86E9 

23 .08E9 

0 

130.87E9 

PHj(dynes/cm») 

G  tdynes/cm*) 

397.14E9 

195.18E9 

468.17E9 

334.31E9 

221.44E9 

50.  E9 

150. E9 

50.  E9 

45.  E9 

150. E9 

G°(dynes/cm*) 

120. E9 

150. E9 

120. E9 

110. E9 

150. E9 

F  (ergs/gm) 

A  (dynes/cm*) 

37.  E9 

10.  E9 

37.  E9 

25.  E9 

10. E9 

5.E7 

0.2E9 

5.E7 

5.E7 

0.2E9 

B 

0.75 

1.0 

0.75 

0.75 

1.0 

Ymax(dynes/c“4) 

0 

3.E9 

0.5 

3.E9 

0.5 

3.E9 

0.5 

3.E9 

0.5 

3.E9 

0.5 

r~i 

0.4 

0.4 

0.4 

0.4 

0.4 

^ambient, 

(ergs7gm) 

1.E8 

1.E8 

1.E8 

1.E8 

1.E8 
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